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IT. INTRODUCTION 


A. BACKGROUND 


Interests in digital image processing methods comes fron 
two principal areas: lmprovement of pictorial information 
for human interpretation, and processing of sc2ne data for 
autonomous machine perceptions. 

Since the second decai2 of this century the interest for 
image processing has becom2 strong, although the theory was 
available at that tine, the absence of processing power 
(speed and bulk storage) delayed the development. Only in 
the earliers 50's, with the revolution of the semiconductor 
industry which speeded up the computers and at the same time 
reduced drastically the storage media cost, the image 
processing technigues flourished. 

Today fast devices and virtual memory comouters allow 
more efficient and simple implementation of algorithms 
without the "complex" and nonuser friendly enviconment like 
overlaying and segmentation of the image. 

The Figure 1.1 shows a typical image processing system. 
The fundamental operations of such a system fali into four 
Major categories: digitization, processing, storage, and 


Grsplay. 


Be DIGITAL IMAGE REPRESENT ATION 


The image signals can be represented intwo distinct 
Groups. The color image and the monochromatic image which 
will be referred throughout this work. 

Any generic image (monochromatic) can be redresented by 
a function £(x,y), where x and y denote spatial coordinates 
and the value of £ at any point (x,y) is proportional to the 


brightness (or gray level) of the image at that doint. 
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Typical Image Processing Systen. 
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Roms cemas Gat al Dpgeescunma sand Storage “a2e concerned 
we need to have a digital representation of images. The 
image £(x,y) which has been discretized in spatial coordi- 
hates referrel to as image samples, is guantized as gray 
levels. Therefore, a digital ilmage can be redresented as 
Medes at grii (that may be uniform or nonuniforn). 

The rows and columns injices locates the quantized gray 
fever at that point. Each element in this griji is calleda 
picture element or simply pixel. 

Although the number of rows, columns, and gray levels is 
feommerary, lt 1S a good policy to work with power of oy due 
to the simple reason that we are processing images in a 
binary machine. 

Throughout this work, a standard (256x256) square grid 
will be used to represent even saapled images with the gray 


level quantized in the range from zecto to 255. 


C. IMAGE PROCESSING TECHNIQUES 


il tEOduc tLOR 


Basically, the technigues HOE arg) ta L image 


processing can be classified into four main categories: 


Image digitization 
Image enhacement 
Image encoding 


Image segmentation and representation 


The above technigues can be implemented in two 
different and highly competitive approaches: the space 
domain approach and the transform domain approaca. Both have 
advantages and disadvantages, and if we consiier the image 
processing as a whole, there is nano clear bouniary between 


each other. 
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For a particular applicacton,) 23 eye O ae o eee 
specific work we have to J0, and “tae Set Of aica “weeeage 
available, we can select the approach among the Space or 


transform domain techniques [Ref. 1]. 


omai 


2- Spatial Domain Tech nmeues 


Spatial domain techniques are based on gray level 
mapping. The mapping used depends on the kind »9f the image 
processing that we are willing to do. sonzales, [Ref. 1] 
describes these mapping technigues in a comprehensive way. 
Histogram equalization, contrast stretching, Smoothing, 
sSharpenig, filtering, contour encoding are classical Space 
domain techniques. More recently [Ref. 2] there are some new 
ones, like image erosion and image dilation, sv«hich can be 


combined to generate opening and closing eff2cts in the 


lmage. 

The mathematical basis for these techniques are 
algebraic manipulations such as vector aeniG, matrix 
operations. 


3. Transform Domain Technigues 


Among the transfora domain represenzations the 
important ones are Fourier transfor, Walsh transforn, 
Hadamard transform, singular value jJecompositiosa (Hotelling 
transform), Kornhuen Loeve transform, etc. [{Ref. 3]. These 
transforms use different discrete basis functions such as 
Fourier basis, cosine basis, Hadamard basis, to reveal 
certain features of the signal. The assumptions which 
Support these methodologies are that the system is linear, 
Shift invariant, and that the signal is uniforaly sampled. 
This last constraint, as w@2 will see in this text, is quite 
Simon 9 drawback for the transform methods in some 


applications. 
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pee -)eSe transkorm domain. rerpseasent ations. are 
Gecived from the basic space domain uniformly sampled gray 
Peameis OF tio Signal, the accuracy of these transforms 
domain representations are only as good as that of uniforn 
sampled gray level itself. Therefore, like the traditional 
Space domain representation, the transform donain ones can 
Beem oe ~Viewsd aS a inear approximation of the original 


signal by using different basis functions. 


4, Nonconventional Technigues 


In order to overcome some difficulties 9f the tradi- 
tional approches for image processing, Such as the uniforno 
samnpling constraint, it is desirable to seek an alternative 
Signai representation from which employs a conplex basis 
functions for better efficiency and achieves a higher 
compression ratio Ehieoug 1 the notion o£ fimict 1 On 
approximations. 

The idea to use nonuniform sampling of image signal 
1s based on tne fact that in overail image pixels have guite 
localized correlation, but over a saall region the correla- 
tion is usually high. Therefore, an adaptive saapling tech- 
Rigue which allows us to have more sampling points in busy 
areas and less ones in doll areas seems to be interesting. 
An approach like this allow us to avoid oversamoling in low 
frequency rejions and undersampling in high frequency 
regions, therefore, reducing the number of necessary pixels 
to represent the image (that is, achieve data compression) 
with adequate accuracy. 

As we will see later that the polynomial spline, 
particularly the B-spline functions, are chosen for such 
alternative representation. 

These B-spline basis have a local behavior that 
seems to match the local sorrelation of the imaj2 pixels. AS 


we will see, there is a tradeoff between the nathematical 


2 


complexity of this New ina9 @seepe seq tee ee eee Loma 


3 
and its efficient and powarful localized bSenavior. 


D. PROBLEM DESCRIPTION 
1. Motivations 


Since the image signal reguires a huge sulk storage 
capability and a wide bandwidth during communications, the 
encoding techniques are a natural way to minimize this 
problen. 

In order to illustrate the problem sone examples 
will be given: 

1) A typical aero-photograohic image £rane | 2567755) 
over square grid will be dealt with throughout this work. It 
requires about 64k of memory, that is, the total memory 
capability of the majority of the microprocessors available 
today. 

ii) The transmission of one of the above frame, a 
line with 1290 baud rate for instance, reguires about 50 
seconds. That is extremely long for nost applications. 

iii) The LANDSAT images are stored on majnetic tapes 
{Ref. 1] where the number of bits required t5 store one 
frame is {234)) (7) (4y = 21 1000000. fhe LANTSA, comlects 54 
frames per day or approximately 2.2 1047 bit p2r year. It 
the data (number of bits) could be reduced by a factor of 
two, the number of data tapes (and the number of warehouses 


required to store them) could be reduced by the same factor. 


2. Objectives 


SS oo eee 


Based on the rationale mentioned in tie previous 
sections, the objective of this work is to explore the 
advantages of image encoding, more specifically data 
compression, by using a nonorthodsx technique. The problen 


of uneven sampling representation will be pursuei. 
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the mathenatical framework is based on the voiyno- 
mial splines functions, more specifically th> B-splines 
basis, which have shown racently [Ref. 4] to be quite suit- 
able for digital image cepresentation, restoration, and 
meter polation. 


3. Formal Specification of the Problen 


Given a digital image represented by its gray 
levels, we are asked to code this image in contours. Fron 
these contours, we select uneven knots based on a criterion 
that gives maximuna smoothness under the constraiat of close- 
ness. These reguirements should take into account’ the 
three-dimensional effects of the image. 

In other words, it may be stated that tae main goal 
1s that the number of selected knots should be Less so that 
ahigh data compression ratio is achieved without intro- 


ducing large errors in the reconstructed image. 


In the next chapter the concept and definition of 


Knot and erroc measurements will be introduced. 


Va 


II. SPLINE FUNCTIONS AND B-SPLINES 


As  ENTRODUCTION 


The objective of this chapter 15 to present an overview 
of the spline functions theory that is the bisis of the 
development of the next chapters. This technigie was first 
formally inte educed as a “Hathenatical ~“cimemwon by 
Schoenberg in 1946 and it has become an important tool in 
various applications, such as, approxination, interpolation, 
and curve-fitting to a set of discrete data points. 

Some of the reasons for the great populacity of the 
spline functions are because they are simple and yet have 
many desirabl= characteristics in the above mentioned apfpli- 
cations. Particularly, the representation of a piecewise 
spline function through a linear combination of normalized 
B-spline basis has some interesting characteristics which 
relates very well to the image model that was seen in 
Chapter I. 


Be SPLINE DEFINITION 


The name spline derives from a draftman's device. A 
spline 1s a flexible strip which can be bended by weights, 
so that it passes through each given point. [t also goes 
smoothly from each interval to the next, accocding to the 
law of beam flexure. The mathematical procedure presented 
here is an extension of this idea. 

A spline function is composed of piecewis2 polynomials 
Satisfying certain continuity properties at the end points. 
AS we can see this is a loose definition. There are several 
formal definitions for the spline function, (Ref. 3], 
{fRef. 6], [Ref. 10], [Ref. 16], although the basic idea of 
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Eee eneormal lefinition reaains the same. Sinc2 we will be 
HeemgG@easmodifying version 9f Dierckx's algoritha in Chapter 
feywhiswnotation is adopted. "Thereforeyythe funstion spline 
PoeeEOrnallysdefamed as follow: a fumction S(x), defimed in 
Rec range asx = b, is callecdsa splinesfunction of degree K 
Meets tCi,eel = ,2,...,m if the following soaditions are 
satisfied: , 


(1) In each interval [t. , «t. Je Ll oye (ay a, 


-1 


b= 4, S(x) 1s given by a polynomial of iegree X or 
less. 
(aay S (x) augers Geb iVvVatives Of “Order 1,2; ...,%=-1 are 


continuous everywhere in the interval [a,b}. 


The order of the spline function is defined as 
m= k*+t 1 C200) 


where mis the order of the spline and k is the degree 
of the spline. 

Example: For a spline of degree K = 3 (cubic spline) the 
order of the spline is m = 4; therefore, we can see that the 
degree of the spline is the high power of tha polynomial, 
and the order of the spline is the number of coefficients in 
the complete polynomial spline. 

The sequence ao 1s called a knot sequence (also known 
as joint sequence or break point seguence for dur applica- 
Gaon here). 

There is a generalizei definition of a spline function 
not only that makes use of the maximum smoothness allowed, 
but it uses plecewise polynomials with less then K-1 contin- 
uous derivatives, which ace known as deficient solines. 

By uSing this generalized definition, the break point 
sequence (or joint seguence) will not be the same as the 


knot sequence anymore. 


1g 


A break point sequence iS a SegueneS Joes item: 
increasing real numbers t, < t,...< t, whereby each break 
point t, is the common joint between two adjacent polynomial 
pieces. 

The knot sequence is a monotonic increasing sequence of 
real numbers XS X,See-5 xX) (M2 Qn), Such that one or more 
knots (x; "S) may be collapsed intS a single Dbreak point 
(ti) 

The advantage of this extended definition 15 that it can 
represent better abrupt transitions and even jimps in the 
function. The less continuous is the function in a given 
break point nore knots will collapse into that point as 
Showh in Figure 2.1. 

As we will see in Chapter III, we are iaterested in 
generating a smooth curve throuyh a set of contour data 
points; therefore, we will use the restrict i2finition so 
that the break point seguence coincides with the knot 


seguence, and they will be interchangeable. 


C. CUBIC SPLINE FONCTION 


According to the definition adopted here, tie degree of 
the spline function is an arbitrary integer number kK. 
However, a plecewise low degree polynomial is usually nore 
Suitable to form a curve through a series of points. The 
advantages in using low-degree polynomials ace that it 
reduces the computational requirements and avoids the numer- 
ical instabilities that arises with high degree polynomials. 
These instabilities can cause undesirable wiggl2s when many 
points must be joined in a common curve. On the other hand, 
however, the low order spline such as the first and second 
order splines, i.e., the pulse function and the piecewise 
linear function respectively, are not smooth. For most of 


the applications they are not suitable either. Tiaerefore, we 
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have a tradeoff between low and high order spline functions, 
based on the above considerations. The cubic spline (degree 
K = 3 and order m = 4) ais advantageous sinc? it is the 
lowest degree space curve which allows a point of inflection 
and has the ability to twist through space. Therefore, the 
Gioee Soline Eunction is a natural choice, andi it will be 
used from now on. 

mie cGubac Spline function S(t) defined ov2r a set of 
mrots (Or break points), (tih , LS represented in each 


interval (t,,t;,,) by a cubic polynomial function. 


S(t): y =a, (t-t;,)3+b; (t-t;) 2+c; (t-t;) +i (2.2) 


more. t, < & < ¢t. : 
1 1+1 


Since it fits at the two ends of the interval, 


Y; = a, (t,-t,) 3+b, (t, -t,) etc; (& Seer ean 232) 


1 1 


pe fi) Pb, (iy tp) eres Chien ti) t9i (2.4) 


Seen, h2tc; h; +d; 


In order to examine the slope and the curvature of the 
joined polynomials, we need the first and the sezond deriva- 


tives. So by differentiating equation 2.2 we get: 


Za 


vy’ = 3a. {(t-t) ) 22h eo eer (25 


yi! 6a; (t-t, teu (2s) 
Now, if we define the second derivative y"* at (t;,y¥;) 
and (Ct: 163i.) Sey =o 2a respectively, we can 
express three coefficients of the polynomial in terms of the 
Given points (t;,Y;) and sthessecond derivativeme of the 


curve at each point. 


(ore) 


oS 

Ci = ¥j-) - ¥y > “2a Seis 
AG 6 

bi = Sj 72 


qi = (3541 see 


Now, from the condition that the slopes 9€ two cubics 


WiuGnesnOd 0 cue (t. ,Y. ) are the same, from equation 2.5 at the 


points t. and t., we have 
yi = 3a;_; (€; ~tji-7 2420, Set, = tS ee (222)) 


AZ, 


eee MaenOms e2.3, 2.9 abd oluging the ssefficients 


Seeessed in equation 2.7 we have the general teco 


eo et cee sone eho Ss, = (2.10) 





The equation 2.10 applies for all given poiats i except 
the boundaries i = 0 and i =n. Therefore, if we write equa- 
move 2s (OR iNneMatrixeform for i = 1,2,...,h-1 wieire S; are 


the unknown we have 


(27a) ae) 


Rei, +b, ) hi, 0 *S by 


h> 2 (hy + h3 ) b3 > De 


| 0 hy-2 2 (h,-2 + Deen Sh 


— 
e 


Dy 


The above matrix has n-2 equations and n unknowns. The 
two extra eguations are given from the end conditions which 
may be given in different form depending upon the specific 


problem, aS we will see. 
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Furthermore, the above matrix DdJecsomes tCLlidiajosa ea ee 
its properties are well known; it can be solved [ast and can 
be stored economically. 

The end sontitions characterize tne different types of 
cubic spline. The common cases are listed below. 


(i) Natural spline 


Si" (en) = 1S Ueno (Pen 


In spite 9f its positive sound name, the natural spline 
Gue to its flatter effects at ends, is not recoumended from 
an approximation theoretic point of view [Ref. 1) ]. 


(Gr ea sp ee 


Si (t) = K, and S*, (t) = & (2413 


where K and K are constant andthe first and last 
Span are covered by a parabola. 
(122) Cyelie or pemtod ie apie 
ie 


S$) (ty ed!) (hy 9 fj = 0,1,2 (2. 14) 


This is quite appropriate to represent closed curves as 
we will see in Chapter III. 


(iv) Antizcyclic spline 


(j) a) 


S Ca) (2215 


1 
— 
_ 
ct 
a3 
eal 
J 
i) 
© 
= 
bo 


ot dp) ; 
S Gn) = oS Gey = 
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(v) Canti-levered spline 


ee) se) tt, ) = 0 (Zeno) 





7142 





Figure 2.1 Multiple Knot Effect. 
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D. B-SPLINE FUNCTIONS 
1. Motivations 


So fac, we have defined and characterizei the piece- 
wise pclynomial spline as a smooth function that fits at the 
break points satisfying some ending point recquir2ments. 

There are some drawback associated with those 
constraints: 

(i) If we are processing experimental iata points, 
most of the time we are not interested to fit a curve 
exactly through this set of points, for the simple reason 
that there is an intrinsic error associated with those meas- 
urements. We rather prefer to fit a curve that presents a 
tradeoff between closeness of fit and smoothness. 

(11) If the number of data points is large, although 
the tridiagonal system associated with the spline function 
can be solved by efficient algorithms [Ref. 7], it is still 
time consuminj and storate expensive. The imag: processing 
work is a clear example of this limitation. 

(iii) The lack of ability to manipulat2 the curve 
fitting, so that we can not alter local behavior without 
altering the whole curve, is a bad characteristic of many 
curve fitting algorithms. 

The concept o£ PLunerion spline reraresentation 
through a linear combination of B-spline basis was intro- 
duced by Schoenberg in 1967. 

The mathematical framework of the B-spline basis is 
the linear space theory. Here we will deal with the appli- 
cations and properties of these basis. For more information 
about these fundamentals see [Ref. 10]. 

One of the first and probably the most important 
application of B-spline for image processing was carried out 
by Andrew, ‘Ref. 4]. In Andrew's work, th2 image was 


sectioned into subimages, each one with a different uniform 
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meet dencttiyeerelatedwetos thewixels distributidmeover the 
Gage, and a Bi-cubic B-spline fiunctionebased on these knots 
was used to represent the image. 


In Chapter III an algorithm based on S-spline to do 


automatic knots selection fron the image contours is 
presented. 
2. B-Spline Definition 


AS we have seen before, we are conceried with the 
Pam y OL SplineS with knots of the multiplicity one. This 
means that the set of given break points coinciles with the 
set of knots. In other words, we are interested in the 
Teaenun shoOtu Spline £unction, i.e., that with continuity 
up to the K-1 derivative, where K is the order of the 
spline. We will be using along this work K = 3, i.e., cubic 
spline. 

For this class of functions, we define the B-spline 
representation as follow: | 

Given a set of knots eae, andes (et Weeeaadi tonal 
knots satisfying 


oy ae eae C= C23} 


nel = th.2 S.-- 8 CR 


Paeewoach 2S further arbitrary. 


Any spline function S(x) of degree K with continuous 


derivative up to K-1, can be represented by 


eg = 2 eM, . 21) (x) (213) 
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where Me) is the normalized B-soline fasis 
and (cain is the coefficient sequence which unijguely repre- 
sents the function S(x) with respect to Whe 3-5) miemeac:s 
(Moga, (80) - 

The B-spline function is constructed ffom certain 
linear combination of truncated power functions, so that 
each basis vanishes outside a small interval. 

Among several definitions for the B-spline (see 
PReL saree as we are going to select the one based on the 
divided difference of the truncated power function. 

Given a sequence of knots Ce), the i B-spline 


basis of order m (i.e., degree K) is defined by 


k +] 
Niger, 1H) =F COG gay TEL) De (Ep eee tiger DY (bed) (2. 19) 
where 
k k 
Gg, (tae) = (tx), (Geek) eee > x 
0 PE: t <ee 
and 


Seats ) f(t) 


th 
stands for the divided difference SE the fEunc— 


Maal - Note that the 
1+1 1+2 


divided difference of (t-x) 1S obtained by £ixang x and 


tion (t) atveae pointes Z, 2 ene 
considering (t-x)* as a function of t alone. 

The resulting number depends, of course, on the 
particular value of x we choose, i.e., the resulting number 
varies as we vary x and so we obtain eventually the function 
(xX )e 


In orjier to simplify the notation, we will use N° 
1 


N. 
lek eet 


instead of Me, (x) as long as the K and the t can be 
1 , = , 


inferred from the context. 
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The normalized B-spline is defined as 


eee) 7 (eo tN Cava 


3. Properties of B-Spline Functions 


The B-spline basis enjoy very nice prop2rties which 
make them guite superior compared with another proposed 
basis. AS we Will see later on, some of these properties are 
very useful for image processing purposes. 

(1) Local basis 

It follows that only K+1 B-splines have no zero 


value at any particular interval [t.. E 


ei Vee ene 


Ae. (x) = 9 Pe kG Ol. xX 2 ieee | (24.21) 
This property suits very well the nonglobal charac- 
teristics of the pixels in a image plane; besid2s, at most 
k+1 (where K = 3 for cubic spline), basis are required to 


evaluate the function $(x) in any point between two generic 


J 


nots x Res 
‘ ! if ell 


S(x) = 2 cimiaer (x) (ee) 


l=j- 


(11) Non-negative basis 


Sa SS 


The B-spline basis is non-negative, 1l.e., 


DE. | (om? ofor t. <x < ti 441 (23-23) 


Zo 


Again this property tS )nteel) SS 6ciat= ie) eee 
processing application once the Signals  vepessen rod maar 
usually light intensity and are always positive Juantities. 

(i111) Partition of the unity 

The summation of all normalized B-spliie basis is 
ecual to the unie, 12e-, 


ee) al (2.24) 


As we will see in the next subsection, the B-spline 
can be evaluated from a recursive relation, egaation 2.25, 
which allows the numeric differentiation and the integration 
to be carried out efficiently by just differing »r Summation 
of the ccefficients respectively. 


For the derivative of a spline we have 


n n 
(j+1) 
DCLG Mi yer (t) J -_ Ci NG ke 4) (2-25) 
where 
ia (k+ 1-35) oe - ei £OE} > a0 
Cee t 
€ for j = Q 


For the integration of a spline we have 


n 
= 
Zea Qa = z a en (t) (2. 26) 
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where 


or = 2 C | (Xe 1+ j me eet wen gt a eee Tl 
i“ k + 1 


DIEtSrESntiation and integration play a1 important 
role in the image processing applications, Sich as edge 
extraction through gradient operator, removal of blur caused 


by linear motion and smoothing. 


4. Recurrence Relation for B-Spline 


———ae SS EE ee 


Although the definition for the B-splin2 basis as a 
Gevyaded diitercence according to equation 2.19 is nathemati- 
Cally perfect, 1t may cause problems like loss of value 
Significance during the numerical computation of various 
Mrererence coefficients. 

A humerical stable formula can be )btained by 
applying the Leibnitz‘'s formula for the cn dnvivaecucdd £fLer— 
Seemomed Product to the particular product. | 


(t-xyk = (t-x) (t-xy**! (Qao7h 


The development of the equation 2.19 can be seen in 


iret. 10 ]. Ee final result Pils the collowinjy recurrence 
relation 
Nef) 7 ees Noe i eee Neel ®) (2. 28) 
- 7c. 
cst i+k-1l 141 
and 
N (x) 1 Ge es) eee 
ae 1 1+1 
0 otherwise 


om 


in the following shapter, tha s-spline bteoms vic 


used as fundamental entities for a data compression model. 


eZ 


III. APPLICATION OF B 


A. INTRODUCTION 


AS was m@ntioned in Chapter I, our objezstive is to 
achieve data compression. 

In the first step, the image contours are geierated fron 
the data over rectangular and uniform grid repr2senting the 
1mage. The total number of contour lines is a function of 
the infomration contents of the image signal ani the numter 
of contour levels chosen (layers). This last parameter is 
under the user's control, and we will see that ia the sensi- 
tivity analysis performed luring the testing phase. 

The rationale to use tontours in the representation of 
image signal is that the contours representation is a 
powerful technique to map three- dimensional data. As 
discussed in Thapter I, a monochrome images can be repre- 
sented as function of two variables and seen as a surface in 
the space. 

Contouring techniques are well known in the past, and 
there are several technigues to generate contoucs. Another 
important point about contour representation is that it isa 
guite efficient way to eliminate redundant information. 
Therefore, it is a suitable tool when we want to achieve 
image data compression. 

In the second step, from the image contours jenerated in 
the first step some "important" knots are selected so that 
the image is represented as data over a set of uneven knots. 
it contains much less elements than the original data set, 
but preserves the important features of the original signal. 
In another words, our image representation with this new set 
reguires much less bulk storage than the original data over 


uniformly sampled grid. 
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The algorithm to select the “important ™ «i0ts wa vonage aen 
Contour is th2 main subject of this chapter. The mathemat- 
ical tools introduced in chapter II (spline aad B=spline) 
are the basis for the algorithm developed here. 

The requirements for this algorithm are stated as 
follow: the selection of the knots should be sich that the 
image could be reconstructed, and that the ercor from the 
original data should be minimized. 

When we jieal with images represented in contours, we 
have two degrees of freedon: the first one is to represent 
each contour (plane curve), and the second one 1S to repre- 
sent the three-dimensional effect, Dae tie relation 
between contours at adjacent layers. 

In the next section, the particular contoating prestenr 
is studied and adapted to the image processing 2nvironnent. 
The following sections deal with the problem 9df selecting 
"important" Knots son caches On tour. These points are impor- 
tant because they provide contribution to data compression 
in representing the image. Informally, their density is 
proportional to the gradient of gray scale, i.e., in regions 
with abrupt changes in gray scale the density of "important" 
points will be high; the opposite is true for dd31l regions. 
The formal criterion for choosing “rmpertale’ srointce. ea 


discussed in Section C of this chapter. 


B. CONTOUR GENERATION 


As we have seen in Chapter I, a digital inage can be 
represented as a function of two variables. These variables 
are the spatial coordinates and the value of the function of 
each coordinate pair is the gray level (quantized) at that 
pol neers 

A contour algorithm determines the seguence df points of 


the image function ina plane, which may be used to draw 
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contour lines representing equal elevation of the surface. 
ghese plane curves through which the image surface is cut 
perpendicularly to the gray level axis are called contour 
levels or simply layers. 

In digital image the gray level is quantized, Lee, 
countable. The criterion to select the number of contour 
levels and the space between them is in some sense arbi- 
trary. In this work, as we are looking for a mor2 robust and 
general approach, the selection method of a number of layers 
and how they are located is not a function of th2 particular 
set of data processed. One of the subjects of this work is 
to study how sensitive is the compression ratio, the quality 
of the reconstructed image, and the computing time with 
respect to the parameters of the contouring proc2ss. 

Contour generating algorithm may be constructed basi- 
cally according to two different methods: 

(1) The tracing method, where the contours ire followed 
from the sam2 starting point until they eith2r close or 
intersect a boundary. The advantage of this method is that 
labeling is celatively easy and that the pen 9f a incre- 
Mental plotter (when we are drawing the contour) does not 
move about randomly. 


(ii) The other classical metho 


(fu 


is performei by dividing 





the image array in small subimages called cells. The size 


WN) 


of the cell i5 arbitrary, and so is the number o£ them. The 
algorithm is carried out by determining all contours inside 
each cell ind2pendently of the other, i.e., each cell can be 
completely processed independent of the rest of the array of 
data. Although in this method the bookkeeping 1s _ more 
complex than in the tracing method, it has sone advantages 
like less memory requirement, and the capability to generate 
contour over a much larger array. 

Perhaps, the most important feature of this scheme is 
its hability to support parallel processing while the 


contour tracing method is inherently a sequential approach. 
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1. Contour Algorithm 3 election 


Based upon a survey of contouring algorithms 
(Ref. 1], [Ref. SJ, “[ Ret.” oF], cur Seequiretens. | arom 
programming environment (see Appendix a), the "“CONTUR™ algo- 
rithm was selected from the NONIMSL Liorary. This algorithn 
is aocontour following type. The algorithm generates a 
contour representationon which one or more (usually more) 
contour levels are extracted, 1.e., given a matcix of numer- 
ical values of z= £(x,y), the CONTUR finds the loci of 
specified discrete values 9f z and saves these l1)dci (contour 


lines) ina file with format specified by the us2r. 


Our task is to extract “Important” Knots (x -amae, 
pairs of spatial coordinates) from the contours. In order to 
do so, we need to have the contours stored in 1 convenient 
way so that our algorithm, which doses the selection, has 
access to these set of data. Other information which are 
hecessary are: 

(1) The density of contours in the neijhborhood of 
each break point gives us indication of how the signal 
energy changes in that neighborhood; therefore, how "“impor- 
tah Gl tse thatekmor. 

(11) The gray value of each contour l2vel (layer) 
and the number of contours that belongs to each Layer. This 
information is helpful in the image reconstruction process. 

From imput data we have: 

(1) The uniform grid representing the image. 

(11) The user's controlled parameters are: 1) "NC" 
that is the namber of levels in which the image is chopped. 
2) "THRES" that is the threshold that eliminates contours 
with less points than "THRES", This last paraneter has a 
snoothing effect that will be helpful when we will deal with 


noisy images. 
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ie oubpuc of the contouring algorithm ace: 

(1) Contour recor3: It is a seguence of pairs (x,y) 
mepresenting the Coordinates of each coutour. [The contours 
are concatenated together and there are two lapels at the 
MegiunangG Of cach one giving the number of points in the 
following contour and its type (opened or closed). 

The Figure 3.1 shows how these sequences are orga- 
nized. An end mark is used at the end of the seguence. 

Pen characte eCr Qatbaxse elt iS a Matrix with the sane 
dimension as the given image containing the projection of 
all contours of the imajye surface. The conto1rs in this 
Meerirx are filly connected with connectivity four defini- 
tion. Each layer is coded by a character startinjy from char- 
Seeer AN (lower contour level) through characte: Z (higher 
coatour level) allowing up to 23 layers. The blaak character 
1s utilized for the background. In this manner, we have the 
Peegeection Of all contours of the image Surface in the floor 
plane. This set of data turns out to be very us2ful when we 
incorporate the three-dimensional effects through the 
HemMsettysOf Contours around each pair of coordinates of the 
contour lines. 

Bae Figure 3.2 gives an example of atypical char- 
eter @md trix. 

(ai Layer record ; lt, vse Simalar t> the contour 
Mapping, as far as structure is concerned. The information 
contained here is about the gray level of each layer and the 
number of contours that belongs to the respective layer. 

The Frgwce 3.3 illustrates 4 typical layer mapping 
with ten layers. 

More information about how these set )f data are 


created, stored, and used are given in Appendix 4. 
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Figure wen Contour Record Format. 
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=< END OF THE RECORD INDICATE 





Figure 3.3 Layer Record Format. 
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Gerona rl SELECTION 


1. Introduction 


Given the contouc information of a1 image as 
discussed in the previous section, we are ready to present 
ame algorithm to fit a smooth spline fwncetion to each 
eemeouL., The Ldea behind this algorithm [Ref-. 10], ({ Ref. 11] 
is to define a measurement for the smoothness of the approx- 
imating splines and another measurement for their closeness 
of fitting so that these two properties, usually contradic- 
tory, could be controlled by a single parameter S. 

The greatest advantage of this method compared with 
others [Ref. 12] is that it allows us to automat2 the deter- 
Mination of the knots. AS it will be shown later on, the 
humber of selected knots iS much lower than the initial 
SemeeouL polnts; therefore, high compression cate can be 
achieved. 


For the sake of better understanding, the algorithn 


is presented in two parts. In the first part, the develop- 
ment is restcicted to a single curve (contouc) without 
considering the effects of its neighbors. In the second 
Dact, that neighboring (three-dimensional) effect is 


included as an extension of the former algorithno. 


ae Pcoblem formulation 


Given the set of data points (x,,Y,), pn 
ieee eES X- 1S in the range [2,5] and y; 1s the value 
of the function at the point x., and the positive weighting 
values Wis i= 1,2...,m, w2 want to determine a spline func- 
tion S(x) of degree K with knots tie Hee kt ieee ht 2 ss. MK, 


such that we have a tradeoff between the following goals: 
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(1) The prescribed value y, @aneu ies page o-- 
enough. 

(ii) The approximating spline should be smooth 
enough, in the sense that the discontinuity of its en 
derivative is as small as possible. 

Mathematically, this criterion can se stated in 
the following way: 


(1) For a measurement of closeness df fit given 


by 
J) = Zw, (y,-S (x,))? (3.1) 
where 
n 
S(x.) = 7 oy ie aoe (Sige) 
ier 


is the representation of a function through the 
normalized B-spline basis, for n selected knots. 
(11) For the smoothing norm aS a measurement of 


the lack of smoothness given by 
- n 
ne) = aie (223) 
r=] 


where d, represents the discontinuity jump of 
the s'*) i()) "a tantameneeeneees 


a = gk) (t_+0)-s‘*) (t -0) (3.4) 


ic 
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OE 


eG 
= Z ans C.. h345)) 


re 
il 
x 


follow: 
minimize 
ac) (20) 
under the constraint 
(ec) < 8 (3.7) 


where S iS a nonnegative parametec under the 
Msem = s Control called smoothing factor. 


Db. PRoeb lem solution 


As we can see, our task is to mininize a func- 
Ponmmimaer da SONStraint. The Lagrange methed is quite suit- 
able for this kind of problem, and it is used here. Fron 


equations 3.6 and 3.7 we san state immediately tiat 





n-k-1 m 
eee iG) cee G2+ ph 7 waelyY —S (xX. \eleShZ2n5S.) (3.8) 
pek+2 | nn 
k k 
Foom the conditions ums & ae 2 = 0 we get 
dp dz 
)(6) = S - 2? (3.9) 
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pie Z = 20 (3.0108) 


AS we are interested in the nontrivial case, 
i,e. p#0 from equation 3.10 we have z = 0. It can be shown 
(Ref. 13] that with every positive value of p there isa 


Single spline S, (x) which minimizes 


. n-k+1 m 
K(pso,6) = 7 d2+pl Zw (S) (x Ny, ) esa (2a 
r=k+2 
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Let's define the function F (p) such that 
n 


m 
p (Pp) = ZS ey ee (in 12) 


n : 
Lew 
From equation 3.9 and our assumptioa on nontri- 
vial case, we have that 


Ee (p) = S (3 ey) 


Therefore, our problem ends up in solving the 
Simultaneous 2guations 3.11 and 3.12. For the first one, we 


evaluate the B-spline coefficients by milniniziny the £unc- 





tion K({p,0,C) for a fi xedea nimee sii o vides hee FOL the 
second one, we find a positive root of F, (p) = S. 

Leading equation 3.5 into equation 3.11 ama 
applying the necessary condition for mininun, Le = 0, 


it turns out that 


(X.Ne Ge] = (3-14) 
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mn 
- = 1a WY. N (xX) i = ly2y eee et 
4 I <1 Lan 1 


Powell ng cqudgtoone oak 21h Matrix form we have 
6C = (34.15) 


where 


ee a 
Se] me Soaiiis lente 

A Vee Xx .)N x 

S| a 1 a i} lak ( i) 
n 

a a ZY, tlhe k (x: ) 
1=1 . 

ee 

wd re (C ne peaieneug ) 


Due to the properties of the locality of the 
B-splines and the assumption that p > 0, the abdve matrices 
A, B, and G are positive, semidefinite, and banded so that 
the equation 3.15 can be solved easily. For more details see 
Martin [Ref. 14]. 

For the solution of the eguation 3.12, an inter- 
active process used starting from the boundary solutions of 
equation 3.15, which approaches to the least syuare spline 
Pumas P approaches to infinity, and to the Weast square 
polynomial P, (x), as ~ approaches to Zero. 


Mathematically, it can be expressed in this way 


n 


F(~) =Z2. w.(y,-S, (x, ))2 (3. 16) 


1=!1 
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F (0) = 2. w. (yee (3.17) 


- L=1 


It can be proved that F(p) iS ContiniouS ene — 
tonically decreasing, and ‘convex “tor p > 0 Sein ores 
Therefore, there iS a unigue root p of F(p) = S. In order to 
determinate p, any rootfinier method can be employed. 

The following rational interpolatioi method is 
caosen due to its fast convergence and simplicity. 

W2 start with three values of p 


Pp = 0, Po= = (CE CO aGelee ote ke (316) 


and the corresponding values of F(p) (FF ,F,) 
are used to cale@ulate the parameters u,@y, w soto <a om, 


functaon ofthe Lorn 
R(p) ="Up toy 7 pvtew (33.45) 


By setting R(p) = S, a new approxination p* for 
is found. 


rol 


p* = Dat cae cE naman > ube ee BE aaa Es ) (F.-S) (3.20) 


Eo (payee. 
3 


p* = -(, -S) (F “Fl pop AUR 5) (ies) > ae) (oe 
E, =S)CF =F, Eee a 


diate Poa ae 
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the value of 7}, Py, P; are adjusted according 
to the following scheme if 


eo aD > «6p pF iee2i) 


Pox? 5 . P3e-P2 Pre P* 


sd) that the interval a P_) 1S permanently 
reduced, while 


oo. < F 
eee, fe SPAS Ss (3. 22) 


The criterion to stop the interaction is chosen 


as 
Faeries s & (3.23) 
Because we are dealing with a smooth process, it 
does not make sense to have € too small; atypical value 


GOemeeoy | 1S suggested by Dierckx [Ref. 12]. 
c. Parametric splines 


Because our contours are arbitrary od9lane curves 


opened or closed, the condition i  T 1S not necessarily 


1 
fulfilled. Conseguently, an adjustment should be done in 
order for the theory presented so far be applizsable to our 
contouring problem. The idea to associate each pair of break 


points (X, oY | ) with a monotonic increasing pacameter 2 ; 
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j = 1,2...,0 appears, to be quite teasenamte: AS it was 
mentioned in Section A of thisS Chapter, zZ iS the free 
parameter and X ; and y; are the dependent ones. 

Then, in this new parametric form, #e can write 
the representation of the x and y coordinates using the 


B-spline approach. 


(3) = 264 iB ie (z) (3. 24) 
a] 
S. (2) eZ cy Mi iat (z) (3.25) 
where 
x = Stee) 
yo SS CZ) Aine S17 AaSe 


y 
ace the parametric representation of a contour. 
Analogous to the former case, the binimization 
of a parametric curve can be stated as; 
Minimize: 


A 
7 (S(t. 499-3 lt -0) J2ec Sk? (t.+0)-S'") (t. -0) J2 (3.26) 
= er xj y j y j 


2. iC Sie re (3. 27) 


Then, analogously to the Minimization problen 
3.8, we end up with the B-spline coefficients Ce y and es 
of the two splines S,,) (2) and Sy (2), derivei for positive 


values of p, as the least square solution of the systems: 
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Yu, 2 Co. fi £25) =Vu, x; 5 = peed eesee (Se 25)) 


& EZ, Sr a j = Vee See aren 
pi 
n 
V5 2 S18 uh 5) = Vw, Y; a Sip 2 pee (Seo) 
n 
TF i= k ie J = 1,2,---,n 


when p is given as the value of the positive 
moor Of F(p) = S, now with 


Meee w(x -S (z))2+ (y -S (z ))2] (3. 30) 
ro J eo) lif 0 8 


About the parameter z, Adams [Ref. 16] suggests 
two methods, and Dierckx [Ref. 15] gives four chdices. As we 
are looking for a general case, because we dd. not know a 
priori the shape of these contours, we choose the Dierckx's 


first proposal as follows: 


Z, =) (Sis.0 1) 
2j = aay d 5 = AUR 7: «som 
See = 2. 

where 


J 


a V (dx ,) 2 + Ca (36357) 
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anal 


AX; = “eae (an 38 


A> Ye ae 

Another important observation should be made 
about the parametric spline as far as smoothing of contours 
is concerned. AS we have seen, the sontour may be opened or 
closed curves. For opened contoucS the paranetric form 
applies directly, while for the closed ones a snall nodifi- 
cation should be done. We just add the periodic condition: 
gj) (a) = s'J)(b) J = Ipse-,K dS a CONStEAING, — comer 
system formed by equations 3.28 and 3.29. 


3. Three-Dimensional Case 
a, introdwetion 
SOmcticiion we have seen the application of the 


function spline to fit a smooth parametric curve through a 
set of m breaking points by generating a set 9f ono knots, 
which is much lower than ao. 

The number of knots adijed in each iteraction 
depends on two factors: 1) The number of knots added at the 
last time. 2) The reduction in the weightei sum of the 


squared residuals Fh (0) - F (o), which is always posi- 


mls 
tive by the fact that the ee are only added and not 
relocated. 

For our purpose, besides the automatic Knot 
selection, another important feature of locating knots is 
that it tends to accumulate more knots in the cegions where 
the underlying data is difficult to approxinaate than in 


areas where a smooth behavior is present. 
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be Algorithm extension 


ive ter otuat the GeqlonsS Wit Grai1cer f1ltting 
error are granted with more knots, guides us to a reasonable 
approach to introduce the three-dimensional eff2cts. We do 
this by assigning to each break point a weight taat is func- 
Etomumor the density of the cortour in its neighborhood. 

‘Then, for each pair of coordinates (a Oe 

every contour, we estimate the density of contour lines 
weommdemt and assign a weight to that point in tais way. If 
the density is high, the weight will be high and the oppo- 
Site is true when the density is low. 

The Figure 3.4 shows an heuristics scheme to 
estimate the density of contour lines at any point on a 
particular contour. The strategy works in the following way: 
Gemtered in each contour’s point (x,y), we coint in four 
Saun@eondl difections (EE, N, W, 5) the number of occurren- 
cies of contours at different layers, 1.e., Ne ; NY el 
N, and select the highest counting among the four direc- 
tions as the representative of the density of tie contours, 
cee MAX (N. , N,, Ny WN. )- 

Mice sSlaheln eteuradirectioOm 1S a parareter "DISpP™ 
that is under the user's control. Then, for evecy pair, the 


weight is assigned according to the following formula: 


x 


W-: 


x 


= 40+ N; AW (3. 34) 


We : “aetial weight 

Wo : reference weight eS) ae Se 

N; 3: maximum number of contours crossing the 
current area 

Aw: incremental weight 


at 


It is the subject of this work to measure and 
reveal the influency of the above parameters 11 this image 
compression method and the reconstructed image. 

In order to clarify the coneept che tae dena. , 
measurement, let's show an example. In Figur2 3.4 we have 
the superimposition of the contour record, Figure 3.1, and 
the character MatGix, Figwee ss. 2. The charact2cr matrix is 
composed of fully connected contours represented in letters, 
while the contour record sonsists of the contour represented 
by uneven spaced knots marked with a circle. The contour 
density measucement is accomplished by moving tie isometric 
cross along each contour so that it 1s centered in every 
contour record point (circled character 10 Pugurs Sa) ean 
an information. about distt&abution of GOnOgcS ane ae 
specific neighborhood can be valiated by checking four 
directions spaced 90 degrees apart. The length of each arn 
of the isometric cross is set by the parameter DISP under 
the user's control. AS w2 will see in the testing phase, 
DISP equals five is a good choice. 

For the particular location of the cross in the 
referred Figure 3.4 we have the following measurements and 
computations: N = 2 Nees N 1 N. =0 =; 


E N i S 
therefore, N = 3 assuming Wo = 1 and AW = 0.2, the weight 


for the break point (X. ay is W. = a6 

Another parameter that has a direce affect in 
the selection of the knots is the smoothing factor, namely 
S, defined in eguation 3.13. This parameter sorks likea 
tuning element that controls the fitness and the smoothness. 
As we decrease the factor S, the function fits sloser while 
1t becomes less smooth due to the large number of selected 
knots. The opposite is true when we increase S. 

There is no optimum value tor S in a wide sense. 
It will depend on how close or smooth we want to fit, 
conseguently how much data compression and fidelity we want 
to tradeoff. 


SZ 


Meee TOL S85 an, INCkeasing EUn=b1on of the 
initial number of break points no. In our testing phase it 
was realized that the computing time increases faster for 
values of m of the order of hundreds. The following expres- 
Sion for S in terms of mo has demonstrated to be a good 
compromise between computing time and amount oof data 


compression. 
S = Ym (m MOD 100+1) (37. 35)) 


where m MOD 100 means: m modulo 100. 


In the next chapter, the inverse prosess will be 
followed. PavhITThOam i gbild wall be constructed from the 
compressed set of knots and compared with tie original 


image. 
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Pictoric View of the Density Measur2ment Scheme. 
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Figure 3. 
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A. INTRODUCTION 


In Chapter III the development of an algorithm to 
compress the digital image was shown. Now, we ace interested 
in reconstructing the compressed imaje and compare it with 
the original one. In other words, we want to hav2 a Prdelity 
Sei1terion to measure the quality ofthe output. 

From the former chapter we have the following scenario: 
a set of knots unevenly spaced lying along the image 
BentourS. ihe density of knots is higher in. dusy recions 
(where faster changes in gray level occur) then in doll 
areas. In order to reconstruct the image ovar (uniforo 
grid), it reguires an algorithm that can generate a surface 
from a scattered data points. 

A variety of algorithnic ideas aas been i2veloped for 
this kind of problem or closely related problems. A survey 
was done inthis area, j;[Ref. 17], rCl esos (uncer aug], 
[Ref. 20], [Ref. 21], [Ref. 22], and two of the referred 
methods based on different approaches were chosen. 

Before going to the description of the select methods 
and reconstruction process, it is necessary to introduce the 


fidelity criterion measurement. 


Be. RECONSTRUCTION FIDELITY CRITERIA 


In image processing system the output image sometimes is 
the final product to be viawed by people or it is the input 
to another image processing unit. The concern is that the 
human viewing mechanism has guite peculiar characteristics, 
[Ret. 1]. Thus, what is suitable for a person may not be 


(and usually is not) suitable for machine and vice-versa. 


a 


Therefore, it 1S reaSonabl= to define two £idelv., cepecr sa: 
the first on2 is based on the human judgement and it is 
called subjective (qualitative) fidelity criteria. The 
second one is based on mathematical relationshid and it is 
called objective (quantitative) criteria. In 91r choice of 
the objective fidelity criteria, we will look for those 


which resemble the human subjective viewing more. 


1. Objective QatantiCative) hese rtey Gerrerrs 


In a image processing systen, like Figure lav waa 
accepted definition for pixel error is given by the differ- 
ence between an input pixel gray level and the sorresponding 
output pixel gray level. Formally, we say that for a given 
input image £E{(x,¥)) as@iietianed yin ey ehagtere IF the error 


e(x,y) at each pixel is given by 


e(x,y) = f(x,y) — 9{x,y) (4.1) 


where g(x,y) is the output image signal; (x, ¥) Gate 
the spatial cdordinates, (27) = 97 1,...,n-1). 

Based on this elemental definition of pixel error 
e(x,y), several definitions of the error which express the 
overall quality of the image arise, and some 9f£f them are 
listed below: 

(1) Root mean Square error 

It is the most popular error measurenints defini- 
CLOM:, not only in image processing, but in many other 
branches of the engineering science. It is defined as the 
Square root of the squared error averaged ovec the image 
array. 


ee = [<e2>]} 1/2 (4.2) 
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where 


n-l n-l 
ae 7 yee- (x,y) 
2 r20 y=0 
N 

Newers ther rotal number of Knots. 

The averaging sense of this definitioi resembles 
that of the intrinsic human visual systen; 1t tends to 
filter out high frequencies and introduce smootiing effects 
in the image. 


(li) Yaximum deviation 





This definition, as its name indicates, assigns to 
the error the maximum absolute value of the 9ixel error 
difference between the input and the output. Lt 2s quate 
often used by the people in graphical area. For image repre- 
sentation by intensity of gray level, it is not juite appro- 
priated because it can give measurements disimilar to the 
subjective fiielity criteria, due to the filt2cing proper- 


ties of the human eye. 


Semen {L(x -y ) = 9 (emery (4.3) 
MAX 1 J : J 
ia 0 pilipwa«y D> | 
(iii) Relative error [Ref. 4] 
It is ilar to the root ean square error. The 


difference is that the normalization is done not by the 
number of pixels N2 but by the square of the input signal 
12 Or ee 


= 
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= 
t 
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Se. CE (x,y) - g(x,y) J? (4.4) 
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(iv) Root mean square of sig atlo 


£ mean a. z 

If we interpret e(x,y) as the noise iatroduced by 
the system at each pixel, we can see the output image as the 
input signal added with the noise e(x,y). Then tae root gean 


square of the signal to noise ratio is defined by 


(SNR) (4.5) 
RMS 





2. Subjective (Qualitative) Fidelity Criteria 


Contrary to the pr2vious section, we can not expect 
here mathematical criteria by the simplest fact that we are 
not dealing with deterministic matter, pathac Ice 
psychovisual affect. Nevertheless, there are som2 character- 
istics of the human visual system which are general, such 
as, the logarithmic senSitivityvesto light “Gneemamty > t7e 
ability to differentiate about 30 gray levels and about 120 
colors. 

The integration effect 1s another feature of the 
human eye that tends to filter out the abrupt changes; that 
implies the introduction of a smoothing effect. 

One absolute scale was developed by Pan2l 6 of the 
Television Study Organization [Ref. 26] and it i5 reproduced 
here: 

(ij Excel Wem: The image is of extcemelly high 
quality, as gosod as you could desire. 

(a eine: The image is of high quality providing 


enjoyable viewing. Interference is not objectionable. 
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(j21) Magginal: The image is of acceptaple quality. 
Interference 1s not objectionable. 

(lv) Inferior: The image is of poor guality but you 
could watch it. Jbjectionable interference is iefinitively 
present. 

(v) Unusable;: The image is so bad that you could 
mote watch it. 

The measurement should be done among a reasonable 
number of persons and the results should be averaged. 

In the next sections twO image reconstruction 
processes will be presented, namely thin plate spline and 
triangulation approaches. 

In the next chapter, subjective and objactive meas- 
urements will be performed on the output images trecon- 


structed by the mentioned approaches. 


Geeeide Ne PLATE SPLINE METHOD 


Hite ackgqround 


The thin plate spline was first mention2d by Harder 
and Desmaites, (Ref. 23). Franke, (Ref. 20], reported 
encouraging results from an algorithm for smooth interpola- 
tion of scattered data by local thin plate spline. The pres- 
entation of Franke's algorithm and its application for the 
reconstruction of the compressed image produced in the 
previous chapter is the subject of this section. 

The algorithm is based on a weighted sun of locally 
defined thin plate splines, andthe final result is aC! 
surface. Despite of the similarities of the nathematical 
relations and notations described here with thos2 of Chapter 
III, we will follow as close as possible Franke2's notation 
because as we have seen, it is not the uniyue solution for 
the problem, and the image reconstruction phase is function- 


ally isolated from the previous steps. 


ao 


Before introducing Franke's 3.367 os hoes 
Convenient to define some terms and set up some aotation. 

Grid: It iS a composite structure 9f elements 
(nodes) arranged ina two-dimensional fashion, Suca that 
rows and coluans intercept orthogonalwise. A uniform sguare 
grid is a lattice with some number of rows and columns 
equally spaced. 

Lattice: same as grid. 

Node (i,j): It is the point where the i-row and the 
j-column cross each other. 

Rectangle Ry: It is defined for Gasca node sia 
not belonging to the border of the grid. It comprises the 
region inside the polygon defined by the four vectices (i-1, 
3-1), (LEV 1) ee eee 

Subrestangle: It is the region defined by two adja- 


cent rows i and iti crossed by two adjacent columns j and 


Ate se AS we can see, each rectangle R coltamms Lows 


Deg 
Subrectangles. The Figure 4.1 shows a section df a uniforz 


Grid with the rectangieg and a subrectangle highlighted. 


Let's now state ene fornal definition of the 
problem:  gGiven™a set of N@scattered pomnts (x7, 77) ee 
l,-«e,N, a lacally defined £Lunction F(x,y) whieh @ias 
property F(x,,y,} = £,, XK = 1,-2.,N, as )to bemeonstrucmoa: 
and from this function we want to generate a aniform grid 
through the evaluation of the function F(x,y) in the nodes 
of a specified lattice, In Our Particular case fere the grim 
has the number of rows "Nx" equals the number of columns 
"Ny", and it 1s Nx = Ny = 256. ‘Therefore, F(x,y) 1s to be 
evaluated in 65536 nodes of the grid. The parameters of the 
function F(x,y) cepresent the discrete spatial coordinates, 
rows, and columns of the uniform sguare grid respectively. 
In our particular case, x and y are integers zanging from 


zero to 255. 
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Figure 4.1 Typical Rectangle and Subrectangle in a Grid. 


The overall interpolant is of the form 


n n 
F (x,Y) oy 7a Wome (x, Y)Q; j (X,y) (4.6) 
ies : 


De 


where F (X,Y) is the value of the funzstion (gray 


level) at the node (X,Y) defined above. 


The summation is performed over the 
The details of this 


"me@mles Of al 


auxiliary grid (dummy variables i,j)- 
auxiliary lattice will be given in the next section. 
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Wij Ke) Sol we yl, dl = eee are the 
weight functions relative to the grid node (i,j) and evalu- 
ated at the point (x,y) that represent the discrete spatial 
coordinate of the output image (uniform square 256x256 
Geis). 

Qi,j; %&-y): analogous to the weight function, it is 
the local thin plate spline approximation relative to the 
node (i,j) and evaluated at the point (x,y) enjoying the 
property that Q(x,,y,) = £. We will see mor2 on weight 
functions and thin plate spline later on. 

The function F(x,y), equation 4.6, has tae following 
properties: 

(cn) nterrolatio 
Ne bss where f, is t 


{Hi 


ins F (Xk + Yk) = “=, . ne 
he given data point (gray level) 
at the point (x,,y,), and N is the number of giv2n points: 

(11) Smoothness: F(X,y) is at least 1S smooth as 
the Wi (x,y) and Qi, (XY); usually ©! continuity 1s guaran- 
teed for all points. 

(iii) Local dependence on the data: F Om a (fred 
(x,y) the weight function does not vanish for at most four 
terms in the Summations defined in the equation 4.6, 1.€e., 
Wi j(x,y) have a small support. 

Similar to the Bespline discussed in ‘chapter III, 
this property is highly desirable because of itsS implication 
in reducting the processing time during the evaluation of 
F(x,Yy)- 

For now it iS important to bear in mind the scenario 
that we have: a set of scattered data™ (x,y 7fi, ), K== 
1,..e,N, and two uniform square grids constructed over the 
region containing the given points. The auxiliary grid where 
the thin plates are defined and weighted , ani the output 
grid representing the reconstructed image whose nodes (x,y) 
are the points in which the functzron ' (7 server one ee 


is evaluated. 
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2. Selection of the Height Function 


MeGOGaing tO the previous Section, it was stated 
that the weight functions have a local support. Although 
ieee, 15 RO theoretical constraint in the Saape of ° the 
support regions, in practice it is common to use rectangular 
regions, whith simplifies enormously the problem without 
loss of generality. 

The use of support regions which are restangle have 
Specit1c advantages in terms of controlling the. number of 
Support regions in which a particular point (x,y) lies, and 
of course it simplifies th2 determination of those regions. 

From these ideas the notion of auxiliaty uniforn 
square grid comes naturally. It 1S construct2?d over the 
given set of points (Xr, ft). K =9i poe py cl mit S Ory eCc= 


tive is to determine the region (rectangles R; where the 


+) 
thin plates will be evaluated and weighted. a aoe daa ly 
feces Gonstructed in the following way: the us2r supply an 
estimative of the average number of points pec rectangle, 
namely "NPPR", and the algorithm using the number of given 
points N ani NPPR evaluates the number the number of 
internal rows and columns that composes the auxiliary grid, 
n, and ny respectively. 

The equation 4.7 #£4x2Shows the described relation 
between the two input parameters N, NPPR and the refinement 
Ometne Grid, n, = n 


X y° 


V2 
n, = NINT[ (4N/NPPR) -1] (4.7) 
where the function NINT means closest integer. 

The weight function W, j (Xe¥) is locally defined for 
each internal node (i randeitSesippoet 15 tie rectangle 


ae The auxiliary uniform square grid with n,+2 rows and 
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eee columns Can be reénresented by 2 he Sosa - eee 


seguence cf [ows 2nd comunns sages r ee 


(4.8) 
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n+l 


Yo < 7, < @ee#¢ < ye yon 


where each pair (7) represents one node of the 
auxiliary grid. We call this sequence auxiliary Lattice node 
sequence. 

The tilde notation was introduced to differentiate 
the auxiliary grid nodes from the given points, spatial 
coordinates (X,, Y,) of selected Knots and the sdatial coor- 
dinates of the output grid (x,y), the reconstructed image. 

Let now 4, (Ss) = 1-3s2+2s3 be the dernite's cubic 
Satisfying H, (0) Sooty 2 Cee B, (0) = 0, E, (1) = 0, as 
shown in Figure 4.2. 


H(5) 
0.00.8 0.2 0.3 0.4 0.5 0.6 0.7 0.8 O.8 1.0 


0.99.f 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.6 


Figure 4.2 Hermite's Cubic. 
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Based on the Hermite's cubic polynonial and the 
auxiliary lattice node seguence, we define piecewise cubics 
with Ci continuity, which do not vanish on two adjacent 
iieervalis, anid satisfy: 


(4.9) 


ad 
ee! x,< < x 
3 Era , Ase >. 49 
Vi (x) { 0 ms. S< ane 
~ ~J 
pewtee (xh og X5.,S ¥ < Xx; 
Ja eS ce 
xX: ~ 1 aA 
0 1+ Pig Ki No 
S 
Vnx (x) 0 7 x < yt 
eri x) Me Xx SEX, 
1 5 ee 


am) 
~ 


Boece jets didlos then.) we cah finally define the 


weight function as 


Wi ey) = V; (x) 0 | (y) (4.10) 
fiat eo cg Nh 
Bee) Vgieas 7 


As we can see fron the definition of equation 4.9, 


the support of the weight function except for th boundaries 
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of the lattice is R. » The referred equation £senalizeamaae 
r J 
local behavioc of the scheme aS was Dentioned 2:arly in this 


Sect Lon. 


3. Thin Plate Spline Basis 


=== aE 


The theoretical developments of the thin plate 


spline was done by Duchon, (Ref. 24]. Analosjous to the 


polynomial spline, it was developed to fit a surface to 
Gectain points Dy eapoerederon of convenient locaivzed 
weight. 


The thin plate spline basis is defined as follow: 


o(x,y) =2 7A a Logg, +. a0+yb toa) (ireed 1) 
(ce : : 4 
where 
I= {k:Q 1S to take on the value f at Cy 9 
x K 
d‘= (x-x,)@ + Oya) and the coefficients Av, a, ey 


and c are determined by the solution of the following time. 
system. 


(4 sen 


z A, ds logan © at bee 


EI y | ie 


(7 yo (X.Y; ) 


Le i AMG 


ce 
Z. 

Ax. = Q 
KEI ee 
2 Ay.= 0 
can k Yk 
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This system of equations has some PiteE2ating prop- 
eels, It iS SyMmetric but not positive jefinit2. The three 
last relations of equation 4.12 has a geometric effect that 
Suppress all terms that grow faster than linearly as the 
Geeteance OL the interpolation points is increased. 

The Figure 4.3 shows the rectangular regions covered 
by the thin plates. Each subrectangle is covered by at least 
four thin plates and the final value of each point (node of 
the uniform sguare grid) is the average of the four local 
computations. 

The implementation of this algorithm cai be seen in 


Appendix A and the Fortran listing is provided in Appendix 
Cr 
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YU CEpeecOMeLEPELY COVERED 
/A WITH FOUR THIN BEATLES 
Fagure 4.3 Rectangular Regions Covered by the [hin Plates. 
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D. TRIANGLE-BAS ED TNPERPOD SE Pie Nees Ou 


1.” Int roiuegievon 


The trciangle-based interpolation method was chosen 
as an alternative approach to reconstruct the compressed 
image generated in Chapter III. The rationale to select this 
algorithm is that it is conceptually different fom the thin 
plate spline nethod seen in the previous section, and it can 
be used for comparison purposes. In Chapter VY, the testing 
phase, the rasults from both methods will be analyzed and 
compared for identical set of input data. 

The algorithm 1s conceptually simple. It consists of 
two main steps. In the first -step it builds 41 triangular 
covering the convex hull of the given set of (< 77S) seca 
using the (x.,y.) data points as vertices of th2 triangular 
cells. In the second step it evaluates an arbitrary point 
(X,Y), 1n our Case a (256x256) uUhlform Ssquaged g21d, by Spa 
doing a linear interpolation inside the correspondent 
triangle where each point (x,y) lies. | 

As we can anticipate, this method genzrates aC? 
surface on the space. The Figure 4.4 shows the projection of 
a typical set of triangles in the x,y plane, where x and y 
are the spatial coordinates. The triangulatioa algorithn 
used here was implemented by Franke, [Ref. 25], using the 


ideas of Akima, [Ref. 17}, and Lawson [Ref. 18]. 
2. Outline of “he Tera nowlati oomsewnod 


A step by step procedure to construct a iniform grid 


over a set of scattered points (x,y) 1S given by Lawson in a 
comprehensive paper, [Ref. 18]. He transcridt here the 
outline of the proposed method: "The critercLon of ime 


maximum angle is used. According to this crit2rion when a 
set of four points used aS vertices of a guadrilateral with 


each internal angle smaller Phabe (eee ee let he qJiadrilateral 
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Memeonvex. One chooses out of two possible ways of parti- 
meomeng the guadrilateral into a pair of triaagles. The 
Pae@itaonhing that Minimizes the minimum interidr angle of 
the two triangles will be seiected. 

In triangulating the x,y plane, we first connect the 
closest pair of points. We next add a point segientially in 
ascending Order of the distance from the midpoint of the 
SlOsestepalreot points. This ordering in addinzj new points 
assures that 2 new point to be added always lies outside the 
polygon constructed with the old points. The new point lies 
outside the circle that is centered at the midpoint of the 
closest pair of points and passes through the most lately 
added old point while the polygon lies inside the circle. 
Each time anew point is added we construct trinagles by 
connecting the new point with the old points that are 
visible from the new Dent and, whenever necessary, 
Nexchange" triangles. 

The second step of the process consists of linear 
Mereepolatwon Of each point (x,y) o£ the uniform square grid 
output with the corresponient triangle constructed in step 
one. In Appendix A some comments on the computer implementa- 
tion of this algorithm is presented, and in Appendix Ca 


program listing is provided. 
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Ve. TESTS, MEASUREMENTS AND CONCLUSIONS 


A. INTRODUCTION 


ifipechapter I the general concepts of digital image 
representation and processing were presented. I1 Chapter ITI 
an overview of the spline functions and partizularly the 
B-spline representation were done. The mathematical proper- 
ties and the relevant characteristics whicao make the 
B-spline a convenient and a powerfull tool to work within 
image processing were carefully explored there. In Chapter 
III the mathematical franework developed in Chapter II was 
used to solve our problem of generating data compression by 
selecting uneven knots from the contours of a given digital 
image. The highlights of the algorithm used thece were that 
the knots are selected in an automatic fashion, so that a 
trade-off between smoothness and closeness of fitness 15S 
achieved for each contour. The relation among neighboring 
contours, namely three- dimensional effect, is iccounted by 
using a weighting function that is related with the density 
PeaecentoursS 10 the vicinity. of each contour. point. The 
development in Chapter IV was towards the reconstruction of 
the image from a set of scattered points generated in 
Chapter III. Two different approaches, namely triangulation 
and thin plates spline, were used in the recsonstruction 
phase. 

In this chapter now, a series of experiments is carried 
out to validate the whole process developed in the former 
Chapters. A set of four selected scenes covering a broad 
range of featires are tested. A criterion for azceptance of 
the process based on performance, reliability, and robust- 


ness is presented. The sensitivity analysis of the systen, 


et 


with regard to the various user's controlied parimeters, is 
done and the results are checked against ouc acceptance 


CEiltechr Lon 


B. SENSITIVITY ANALYSIS 
1. Qbjective 


The rationale for a study of the sensitivity of the 
whole image processing system with respect to the user's 
controlled parameters is the following. We can, based on 
our acceptanc2 criterion ts be developed, judge how good is 
ouc model. AS we have seen in Chapter I, our goal is to 
achieve data compression within a reasonable processing 
time and reconstruct the image within acceptable error 
compared with the original one. The factors that influenced 
the performance of the image prosessing system are the 
intrinsic characteristics of the algorithm described in the 
former chapters, and the user's controlled paraneters which 


are the subject of our Study here. 


A desirable image processing system shoild have the 
following attributes: 

Generality : the system should be applicable toa 
variety of images or at least to a complete class of images. 

Predicability : given an image to be processed, the 
performance of the system and the characteristic of the 
Output product can be anticipated based on tie previous 
knowledge of the system. 

Robustness : the algorithm which the system is based 
upon should be strongly formed. It should have 2 solid math- 
ematical fram2work over which some heuristic approaches can 


be adopted. 


vez 


Perfotmance ; the processing should be accomplished 
in a reasonable time with the required precision. 

Simplicity : the system should be simple to under- 
stand and simple to operate. 

modularity the technique, divide saad Weonguer, 
should be used to design and implement this digital image 
processing system. A system easy to understand, expand, and 
Maintain can be derived using the modular approach. 

Gracefuli degradation : an image processing systen 
should not collapse with the small anomalies or pseudo 
anomalies in the image to be processed. It should have 
alternatives to patch small perturbations on the system 
automatically. For example, in the process of finding the 
cootours we have two categories: the opened cont) urs and the 
closed contours. Eventually, if we have asaddle shape 
negaen, weecan end up with a contour that 1t as neither 
Seetiemw Or Closed. The algorithm should be able to handle 
this problem in such away that it does not crashes the 
program and will not have unexpected result. Another case is 
in the evaluation of the roots of an equation. If the 
limited number of iterations is reached before the required 
precison is obtained, a compromise solution should be 
achieved instead of occurrence of a catastropiic failure. 
Like these two examples w2 can list nany others. Many of 
these goals are in conflict with one another, this trade-off 
1S a necessary action to be taken. 

Next, the input variables (user's controlled parame- 
ters) will be described followed by the output parameter. 
The sensitivity of the system will be studied based on the 


variation of the input variables. 
3. User's Contrclled Parameter 


in oObader to simplify the sensitivity analysis of the 


image processing system developed in this wock, ceo S 
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advantageous to explore its modular organizatioa and divide 
it into stages. The natural division at the high2st level is 
to have two phases. The first phase concernes the data 
compression and the second deals with the image reconstruc- 
tion from the compressed data set. The first phase can be 
further subdivided into the contouring phase aad the knot 
selection phase. With this structure in mind, w#e can study 
the user's controlled parameters and their influence in each 


phase of the process. 
as Contouring Phase 


In this phase the contours are extracted froma 
given image represented by a uniform (256x256) square grid. 
The parameters and actions under the user's contcol are: 

1) Filtering : alow pass filter i153 applied on 
the raw image for smoothing purpose. As we know, we may have 
an image contaminated with noise, 4nd particularly, thaceas 
the case of a@rialphotoyraphic pictures. The causes of this 
spurious effect that may appear in digital imajy2 are due to 
poor Sampling system or transmission channel. [n_ our case 
here, a movinjy average filter is applied by the aser who can 
Optionally do the filtering operation or bypass it. The 
filtering operation is the first possible treatnent oon the 
Signal (input image). 

1i) Contour layers : the number of sections in 
which the image will be contoured is called "NC". According 
to the constraints listed in Appendix A, this parameter can 
vary in the range from two to 20. In practice, as we will 
see later on in this chapter, a typical value foc NC is ten. 
Small values tend to introduce large error in the recon- 
structed image and large values tend to deg2znerate the 
compression ratio. 

14) Threshold : this parameter is called 


"THREES", It fixes the elimination threshold that scrappes 
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the contours with less points than THRES. THLS parameter 
like the filtering has a smoothness effect. Th2 difference 
Pomtina@t it “dots after the contours are generated while the 
filtering oftion is applied before. Besides its low pass 
EPeeering Characteristic, the "THREES" parametea: can reduce 
the amount of data to be stored. Therefore, it is a positive 
factor if we want to have high compression ratio. On the 
other hand, the low pass filter has side effects; removal of 
the high frequencies in a signal means degradinj the sharp- 
ness of the eiges. Thus, we must trade off these features 
permepetter results achievement. The range fron five to 10 


Pome "TARES" has demonstrated to be good. 
b. Knot Selection Phase 


This phase corresponds to the second step 
towards the data compression, and constitutes the main core 
step of the image processing system here. The algorithm to 
perform the knot selection task was developed in Chapter 
III. Its main objective is to extract the "strategic" knots 
of the contours generated in the first phase. Intuitively, 
the idea iS to choose more knotS in regions where the 
commourse beni with Small radius of curvatuces {planar 
effect) and/or in regions where the density of the contours 
4s high (three-dimensional effect). Mathematically the 
problem is a minimization of a function under constraint, 
j.e., a Lagrange tyre of problen. In this phase there are 
Six parameters under the user's control which ace discussed 
below 

1) Reference weight (#0) : it is assigned to 
each contour point (x,y). The typical value for #0 is one. 

ii) Incremental weight (4W) : it gives each 
contour point the weight according to the density of 
contours crossing the neighborhood of that referced point. 
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iil) Neighborhood span (SPAN) : it defines the 
size and shape of the neighborhood where the contour density 
measurements will be performed. In order to simplify the 
implementation, a squared region centered at the point under 
measurement is chosen. The size of the square LS two times 
the SPAN plus one (see Figure 3.4), and the typicai value 
for the SPAN is five. 


iv) Smoothing Eactor (S) : i1t is tae parameter 
that sets the trade-off between smoothness on closeness of 
the curve fit. The two bound values for S are *S = 0" that 
results in the interpolation polynomial with no data 
compression at all, and"S =© =" which generat2s the least 
Square polynomial with the theoretic maxiaim of data 
compression. The factor S$ Can be adjusted Sone vcuee, 
between the apper and lower bound. It 1S chosen uSing an 
heuristic approach that gives an acceptable fitiess without 
geopardizing the compression ratio and the computing tine.. 
It is known that the factor ""S" ougnt to be aq ircreasag 
function of the number of points (x,y) in the contour being 
worked out. Several authors, Dierckx {Ref. 5], #=Reinsen 
{[Ref. 27], have proposed empirical ways to choos2= the factor 
S. Here in this development we are going to use 1 generating 
function for S given by eguation 3.35. This has iemonstrated 
good fitness results for contours with both small and large 
number of points, while keeping a quite reasonable computing 
time and producing a substantial data compression ratio. 

v) Tolerance (TOL) : this parameter sets the 
tolerance by which the root of the equation 3.23 should be 
evaluated. Once our resolution is the space between two 
adjacents nodes of a uniform sguare grid (256x256), there is 
no advantage in using asmall value for TOL because all 
selected kncts will be round off to the nearest node on the 
grid. Besides, the computing time is an inverse function of 
TOL. A typical value tor TOL gs" 10) — J 
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vl) Maximum nuaber of iteration (MAXLT) : this 
Bamgameter Sets the upper limit for the number of iteration 
imeene Ssobution of equation 3.23. This parameter is also 
Eomaccamwith thesconputing time during this secondwphase.m A 


typical value is MAXIT = 20. 
Cc. Reconstruction Phase 


In Chapter IV two different image resonstruction 
approaches were used, namely triangulation method and thin 
piate splines method. The first one does not require any 
input parameter from the user except the files of input data 
generated in the previous phase, as it is illustrated in 
Figure A.1, whereas the second one reguires an extra paran- 
eter, namely "NPPR" to be supplied by the user. This paran- 
eter isan 2sStimate of the number of points per region 
(defined in chapter IV). The reconmended valie for this 
parameter is in the range from eight to 30. . 

So far we have described ten user's controlled 
parameters and the respective ranges. It is easy to realize 
that for a set of different images combined with a_ set of 
ten input parameters setting we certainly would end up with 
huge number of testing cases. Therefore, befor2 discussing 
the output parameters and the test itself, it is worthwhile 
Oo Narrow the variability 9f the variables. Using the past 
experience ani some reasoning, we can reduce tie degree of 
freedom by fixing some of the described parameters. The 


following parameters were chosen to be kept constant. 


parameter Status/value 
FILTERING yes 

SPAN 2 

5 VM (M MOD 100 + 1) 
OL 0.001 

MAXINT 20 


w/ 


where M4 1S the total umber of | Ppermro ae: 
contour and the other names were defined above. 

Thus, the effective varying paran2ters which 
will be varied in the overall procedures are: NC, THRES, NWO, 
AW, NPPER. 


4. Qutput Measurements 


Once the user's controlled parameters ar2 presented, 
it is time now to introduce the measurements. The sensi- 
tivity of these measurements with regard to th? paramenters 
is to be studied. Those measurements are: CPU time in each 
phase of the process, storage required, compression ratio 
achieved, and qualitative and yuantitative error measurement 
between the ociginal and the reconstructed imajg2. With the 
definitions and assumptions set we are ready t3 start the 


testing that is the subject of the next section. 


C. TESTS AND RESULTS 


Four aerialphotographical images were selected to test 
and verify the image processing moiel. These test sets were 
chosen to cov2r a wide range of cases, i.e., from a very 
doll image to a heavy busy one in order to check the accept- 
ance criterion defined in the previous section. 

For each of these images the whole procelire, 1. Ce, 
contouring, «not selection, and ~reconstructroa detined ma 
the earlier section, were recorded and the results appear in 
Tables I to Vili and any figures 5.1 tem. 

The Figure A.1 shows the diagram of the i1nage process 
System, where the three phases are representei and so the 
interaction between programs and files. As we have mentioned 
in Chapter I, the test images are represented by a unifora 
sgGuare grid with spatial coordinates ranging from zero to 
255, and this same range is used for the gray level 


quantization. 
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Let"S now Summarize the notation used in tie tables of 
the test cases. 
i) First Phase : Contour Generation 
DPiteeo > Identatication of he =data rile under test. 
NC s number of contour levels in whica the input 
datvemwiatiebe divided . 
Dah ES : threshold that defines the contours to be 
purged. 
TCTCONT : total number of contours generatei which meet 
the threshold value. 
POLeAehe. total enumber or GCSOEdInates Pal i eee) 
representing the total number df contours 
(TOTCONT) . 
MINLEVEL: minimum gray level value found in the given 
data, ranging from zero to 255. 
MAXLEVEL: maximum gray level value found i1 the given 
data, ranging from zero to 255. : 
CPUTIME| : total of CPU time inh minutes required to 
accomplish the current phase )f the image 
processing. 
11) Second Phase : Knot Selection 
EXP : experiment number. 
WO 3; reference weight. 
QW : incremental weight. 
N >: humber of selected knots, (x,y) pairs. 
N' ;: effective number of selected knots. 
Cp, : compression erroc achieved (see Appendix B for 
details). 
TCOM: total time required for the compressioa phase 
(contour generation and knot selection). 
iii) Third Phase : Image Reconstruction 
TRIA : nunber of triangles generated. 
Eays 3 root mean square error. 


Fact : relative error. 


uo 


NPPR : eStimate of nuaber of poiots per rectangle. 


The four study cases will be presented next. The data 
for each case comes in the following sequence: 

- two tables containing the parameter setting and meas- 
urements of each phase. 

- One picture of the original image. 

- Three contour pictures showing the contouc map gener- 
ated in the first phase, followed by two contour maps gener- 
ated in the second phase for two different parameters 
setting. 

- Two pictures of the reconstructed image,  Idne for each 


metnod, triangulation and thin plates respectively. 
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PARAMETER SETTING 
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Figure 5.2 ~Ontour Repres. Before Knot Selection - Case I. 
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Figure 5.9 Contour Repres. After Knot Selectioa - Case Il. 
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Figure 5.10 Contour Repres. After Knot Selection - Case Il. 
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Pigure 5.11 Reconst. Image Via Triangulation - Case II. 
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Pigure 5.12 Reconst. Image Via Thin Plates Spline - Case II. 


90 


| 
| 


f 


O¢ 


i es es 
led 

fev oro 
a eee far oe 
[ore re oe Po 
mf refs 
rae 





uw) 

a 
~ 
— 


La 
Cones, Ol rsaorl we [rue de 0 R°O 
or fcc oe [so feo 
ofc [ tno [oo 

vf fcr re ro 

G 
re 
rer 
= 
Poe 


™N 
= 


oo. 
oO 
ON 


= 





wy) 





LGORL 
a 


HO°R 


Cone 
da. 


NSN 
— 
help 


| i ae ae [ 
pews [oo 
ps] ws] | 
cca [rel [ 
aaa 


WIL 
aoe WOOd Ado N 


© 
— 


a 


— 
Oo 
= 


90° Cl 
Sha, 
NOLL 


—, C.J 
22 oe 
Oo, 
< 
_ 
jos 
4 
ad 
Eo 
oa 
ES 
8 
Poel 
ee 
fi) 
C.J 
a= 
es 
oO 
= 











OLLVYANGUS YNOLNOO +ASVHd LSYIS 


Ted FSITH - TII 9sSeyD 3S9_L 
A dTavi 


EE a LSS a ee ee ee ee ee ee ee es ee re eee eee ee LE er cr re cc ee ee ee oe 
eee 





MWELAdO | THADTXVW | TAADINIWN | YLVd LOL LNOOLOL So dMHL Soon alieeh | 4 


oa) 





| 


e——— 


18 
20 
eae 5 


1G ote 


ole 


NerRVeCeU. LE iE TRIAYCPU JE E 


Z10 


oz 
ee 
3 
= 
Oc 
fe 


IMAG 


* 
s 


i 


Second Part 
THIRD PHASE 
tN SPEATES SPLUN 

55 


CONTOUR GENERATIO 


Test Case III 
PTL ae > THRES TOTCONT | TOTPAIR | MINLEVEL| MAXLEVEL | CPpUTIME 
807 


rm is i rm ee i a eg ce em ep cm cr eg i ee ee mm mmm SS SES 


20 






5087 [2 


49 - | ol LOSS SY 72749) Tea GCwaO 10% 
48 pis | 53 [iro] 7.52) 14961 85 Ps 27 | 7 ac 


TABLE VI 
48 


“See. 
walen| as [wo [wool rol | | 


, 
a ro 


Pbk PUP aoe 


1 
4 


fon bor ee Te [e faofof— 


CEU 
RTE 


alia 
ec 


0.8 Od 12069 50 72 On boud A 


KNOTS SELECTION 


Pa 
x 


maf Popool@ [o [ > foprfiesl on fo 


eas | 
a 


[ice fos fn fe 


pee 


SECONDS PHASE 


0.8 OS Ze 20091 30 15 20m. 7/1 


1.0 O.1 711214 
1.0 O.1 | 124 


MIR.DAT 


)6 | a0 | 0.2 fasrio 


la 
fa 
al) 


[Pe fem 
Peo fow 


[ne oe eee ee er rr cam rm erm rr cm a rr cpr rc ror ee cm er ee ee ee ee ee ee ee ee 





o2 





FPiguee 5.78 Original Image - Case III. 





PARAMETER SETTING 
NC 10 
THRES = § ' 


ua 


Ge 
Contour Repres. Before Knot Selection - Case III. 


Figure 5.14 


oes 





4 i 
Ss a = 
d ; yf «@ 
} : / FR THRES = $ 
re’ / = 1, 
IF... > Pran -b 1 Ow =0 


_— eo. 


Figure 5.15 Contour Repres. After Knot Selection - Case III. 





Figure 5.16 Contour Reeres. After Knot Selection - Case III. 
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Figure 5.18 Reconst. Image Via Thin Plates Spline - Case IIl. 
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Figure 5.20 Contour Repres. Before Knot Selection ~ Case IV. 
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Pagwre™S.22 Contour Repres. After Knot Selection - Case IV. 


Case IVY. 





Figure 5.24 Reconst. Image Yia Thio Plates Spline - Case IV. 
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Be “HESOULT ANALYSIS AND COMCLUSIONS 


Bum ole Easenhical pattern will be used to analyze the 


beomerts presented in the previous tables and pictures. Yirst 


all the problem will be tackled as a whole. The analysis 


tty 


O 
of the intermediate results, and the performance of eack one 
or the three phases will be done separetely. 

The data compression ratio generaily dep2nds on the 
contents of the image. The compression ratio decreases when 
the image fecomes busier. In another words, we have a 
Smailer compression ratio for an image with a Large number 
of features. 

In the four study cases presented here the above conclu- 
Sion can be checked easily. The compression 7zatio ranged 
from a value of about 30 in a doll image as in study case 
Il, toa value of about seven ina very busy image as in 
Sseuay case IV. 

Compared with the current benchmarks ji2veloped by 
Gonzales, [Ref. 1], and Anirews, [Ref. 4], these compression 
ratio values are remarkable. In those referances, the 
typical compression ratio achieved was from two to seven and 
from two to ten respectively. 

It is necessary to recall that the coding aljorithnm used 
by the two listed authors is different from what we are 
uSing here. 

With regard to the fidelity of the reconstricted image, 
as it was expected too, the discrepancy between the qual ta- 
tive measuremant and the guantitative ones (root mean square 
and relative 2rror) is greater for the busiest inages. This 
can be seen if we compare the case study II, tne doll one, 
against the case study IV, the busy one. For the same 
setting of parameters we see that the case stuiy IV has a 
Guantitative error of about three times the case study II, 
while the qualitative measurements in both images give about 


the same guality for then. 
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Compared with the results of the above two references, 
the average quantitative 2rror in the process d2ascribed nere 
was higher than those ones. Andrews result shows compression 
error from 0.2% to 6.3% based on their test inages. Our 
result shows compression error fron 2.0% to 11.0% on the 
images used. These errors are based on two different set of 
test images. It was not found a specific reason for being 
so. Althoush if we compare the amount of deterioration of 
the reconstructed images in a gualitative fashion, there is 
not too much difference between our Samples hece and those 
described in the references. 

If we consider the product of compressio1 ratio and 
guantitative 2rror, and use it as a figure of aerit of the 
image processing systen, we will find that the system 
described here and those of the two references mentioned 
above have about the same value. 

The third factor to be considered is the computing time. 
In order to geteethe touametinemeopme themumage processing] 
system described here, we need to add the three time quanti- 
tity appearinjy in the study case tables. The two figures 
refer to the compression phase itself, and tae last one 
refers to the reconstruction phase. From the result it is 
clear that the computing time is high, even thoigh an image 
with (256x256) is considered a large set of data. The 
processing speed of our <soOmputing environment is not’ the 
state-of-the art. It should be convenient to ramember that 
the time shown in the table rated as CPU time iloes include 
the swapping time required to input and output pages into 
the physical memory, andit is well known that large array 
processing generates a huge number of page faults. 

As we Can see from the tables, the reconstruction method 
used (phase three) contributes heavily in the total 
processing tine. We will see later on, that desending upon 


the reconstructed method chosen, the processing time will be 
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fomeee SCUS itive of Not to the amount of data to 2e processei 
(the selected knots, as it was named before). 

Now, that the overall behavior of the system was 
presented, it is necessary to discuss the performance 
details of each phase. 

The first phase, the contour generation ssheme, have 
demonstrated to be avery efficient and general algorithn. 
From previous tests, not shown in the presented tables, it 
Was noticed that variations of the threshold (THRES) param- 
eter has minor effect on the output results when THRES lays 
in the range from five to 10. As a matter of fact, the value 
of THERES equals to five gives a slight improvenent in the 
gualitative and quantitative error measurements between the 
output image and the original one. The compressijdn ratio and 
the computing time have a slight deterioration as it was 
anticipated. In order to reduce the degree of fr2edom of the 
problem, therefore the testing load, we felt that the paran- 
eter THRES could be kept constant without undermining the 
test purposes. Thus, the value of THRES was kept equals to 
five for all experiments but one recorded in the presented 
tables. 

The number of contours "NC" was the effective varying 
parameter in the phase one. A cange from 19 to 15 was 
covered for all study cas2s, anda continuous improvement 
was achieved in the guality of the reconstruct2d image for 
higher values of "NC" with expense of a continudus deterio- 
ration of the compression ratio and the computinjy time. The 
mentioned range for "NC" has demonstrated to be ideal. 

The second phase, the knot selection algorithm, is the 
core of the image processing system developed here. It gives 
the strategy to select the knots so that a high compression 
ratio can be achieved with small loss in the information in 
the result. 
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Comparing the contours plots generated in this second 
phase with those ones produced in the first phase, it 1s 
easy to get the feeling 5n how the algorithm 4 ompmeee Each 
corner of the contour lines represents a s2lected knot 
(output of the second phase), which is clustered in regions 
where the curve bends with a small radius ani/or regions 
with high density of contour (fast change in gray level). 

So far we have discussed the dynamic of the second phase 
through the analysis of the contours generated there. A 
natural guestion to ask is, how does this schen2 affect the 
final product, 1.e77™ the reconstructed mage, =slcanly , / ase 
will inherit the coupled effect of the three phases of the 
process. In order to isolate the effect of the second phase 
in the reconstructed image, it 15 necessary to fix the 
parameter in both, first and third pnases ani analyse the 
final results. This procedure was done many times and the 
results were recorded in the study case tables presented in 
the previous section. 

From this we can extract the effects of tae reference 
weight "WO" and the incremental weight "OW" on the quality 
of the output image. These two parameters allow a fine 
tunning of the fidelity of the reconstructed image. A set of 
value which gives a good trade-off between the error meas- 
urements inthe reconstructed image and the compression 
ratio is WO = 1.0 and AW = 0.1, which proved to ve satisfac- 
tocy for all four study cases. 

In the third phase the two independent approaches, 
namely triangulation method and thin plates spline method, 
have yielded interesting results. 

The triangulation has Jemonstrated to be a v2ry general, 
predictable, robust, and simple method. It covered the four 
different stuiy cases without any constraint. Regarding to 
the performance, it has demonstrated to be guite sensitive 


to the size of the input set (select number of kaoots). For a 
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Gee meron put data set of about 3500 points, the CPU time is 
Similar to that one reguired by the second rezonstruction 
method, the thin plates spline methoi. But 2S the input 
number of data points increases, the computing time 
ancreases very fast. This iS an intrinsic chararsteristic of 
the triangulation approach and nothing can be done about it, 
once a hull of triangles should be construct2d over the 
given points, and the number of triangles 1s of course an 
increasing function of the number of points. The act of 
constructing triangles and working with it becomes very 
expensive in terms of time if we have a lot of them to deal 


Wacth . 
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In Figure 5.25 is shown the tine spent In” DoOGmemecon] 
SEFUCtTION Tethode, triangulation and thin o.lapes spire: 
versus the nunber of input points to be processed. 

For the quantitative error, the triangulation aethod 
always showed a higher value than the thin plates spline 
method, but the difference is not so significative; but for 
the qualitative measurements, a sort of geometric effects 
due to the triangulation were present in the inages recon- 
structed through this method. | | 

The second method for image reconstruction, the thin 
plates spline, did not demonstrate to be gquit2 general as 
the triangulation method, neither quite predictable in the 
sense that it can overshoot, and the opposite happens in the 
triangulation method; and depending upon the sp2cial condi- 
tion on the jistribution of the data (region with lack of 
data), it can generate high overshoot, thus saturating the 
gray level scale. 

Its performance is very good compared with the triangu- 
lation method; see Figure 5.25. As we can conclide from the 
mentioned figure, the thin plates spline methoi is almost 
insensitive td the number of input data points, which is a 
vecy nice characteristic. The quantitative ercor, as we 
said above, is always smaller than that given by triangula- 
COW. Conversely to the triangulation method, the thin 
plates spline method does not have the geometric effect and 
its fidelity is much better than the first method. 

Conversely to the congenital defficiency of the triangu- 
lation method, for the thin plates spline methoi one can get 
rid of its own drawback, strong overshoot in ar2a with poor 
distribution of points, by doing some small changes in the 
program. The specific changes are out of the scope of this 
wock. For more information see Franke's work, [R2f. 19]. 

For what we have seen in this section, the thin plates 


spline method 1s a better option to reconstruct image than 
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@ae t@iangulatwon method in many aspects, and it is poten- 


tially better in an overall sense. 


Es. RECOMMENDATIONS FOR FOTURE STUDY 


The present study was done concerning the fiture use of 
some of its contents for real time application in the data 
compression area. The results achieved in th2 first two 
phases, mainly in the knot selection one, were very encour- 
aging. The scheme proved to be very robust and in general it 
has been demonstrated to jegrade gracefully. 

As we have seen the process was implemented in a sequen- 
tial fashion. A parailel processing, although nore complex 
=O implement, would decrease the processing tine 
considerably. 

Regarding to the third phase, the image reconstruction, 
only two methods among many possible were tried. It would be 
worthwhile to study this subject more deeply. A method 
developed by Lawson, [{Ref. 18], would be a good candidate 
choice. 
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COMPUTI 


The algorithms discussed in this section ace based on 
the mathematical developments presented on previdus chapters 
of this text. They consist of a set of computer orograms and 
files which are used aS communication intecface among 
programs. 

In the next two sections the structure, constraints, and 
usage of the programs and files are presented in details. 
These programs were written and/for adapted to the progran 
support environment of the VAX 11/780 and VAX 11/7590 
machines under the operating systems V4¥S 3.4. Tie languages 
used were Pascal version 2.1 and Fortran version 3.3. 

The Figure A.1 shows the interaction betw2en programs 
and files in a chronological fashion, i.e., the flow indi- 
cated by the arrows represants the evolution of the process. 
The basic elements of such a block diagram ace programs 
represented by rectangles and files represented by hexagons. 
The ovals mark the beginning and end of the process. Inside 
each block it was assigned an arbitrary name for programs 
and files respectively by which they will be called from now 
on. 

The process aS it was mentioned before consists of two 
main distinct phases. In the phase one our conc2rn is to do 
data compression; from a given set of uniform sampled and 
guantized data representing an image frame, the sontours are 
extracted and from these contours an algorithm based on the 
B-spline properties and Lagrange's minimization criteria is 
applied so that the number of knots necessary Eo represent 
the image is irastically reduced, therefore achieving data 


compression. In the second phase the invers2? process is 


108 


Saeercawollt, loce., tne reconStruction of the imzge based on 
the set of compressed data is attenpted; this was done 
through two different approaches namely the thin plates 
spline method and the triangulation nethod. The final result 
is cecmpared with the original one uSing gualitative and 
guantitative measurements. 

In the next two sections the programs and files which 


accomplish these tasks are presented. 


Aes Ee ROG RANS 


1. Program Name: WDAT1.PAS 


1.1 Objective 
Given apacked file MYFILE.DAT containinjy the image 
data, this Pascal program unpacks the file and rewrites it 


Byued SUttable form to be read by a Fortran progran. 


tezeeinput Data 
File MYFILE. DAT: packed array containing the image frame 


momeeecimes Cxpanded, 1-¢€., a (512x512) grid of pixels. 


1. 2eGUES Ut Data 
File IMAS1.DAT: unpacked array containing the image 
frame converted to an array of (256x256) pixels suitable to 


pe read by a Fortran progran. 


1.4 Program Structure 
This 1S a Pascal program composei sf a singl2 block with 


no procedures or function calls. 
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1.5 Constraints and Assumptions 

The input file is assumed t9 be a pack2d array of 
(512x512) bytes ranging GEOmagero £06 255 0y 12h theme wer 
stripped off. 


1.6 Common Areas 
Net applicable imgtols case: 


1.7 Error Messages 
No data in the file: 1 the anput f£2le is san enet ree 
EOF (end of file): if the file is smaller than expected. 


1.8 Usage 
Compilation: $PAS WDAT1 
Linking : SLINK WDAT 1 
Running : S$RUN WDAT1 


2. Program Nane: THE8M.FOR 


2.1 Objective 

The image, IMAG1.DATI, consists of an array X£ (256x256) 
pixels ranging from zero to 255. This Forttan progran 
extracts the tontours from this set o£ data accomplishing in 


this way the first step in the data Compression ohase. 


2.2 Input Data 

File IMAG1. DAT: unpacked array of siz2=  (256x25a 
containing the pixels of the given image. 

File FORO10. DAT: set of parameters which controls the 


contouring operation. The parameter "NC" is th2 number of 


Paiweiommhich sthe image is to be chopped and "TARES" is the 
tiresmetwd that prunes the Sontours with the number of points 
belowmmelHRES. 


Pao eoutout Data 

Popece tr OR@UC. DAT: Character matrix wap. It 15 an array of 
Character with Size an 2) OX 906} which repr2sents the 
Beoavyection Of the contours with connectivity four on the 
plane of the spatial coordinates (x,y). The background is 
represented by the blank character and each layer is repre- 
Pommecemey edn aipna Character starting with "A" for the lower 
level. 

File FORO19.DAT: number of contours per eash layer and 
the respective value (gray level) of that layer. 


EieboenOno20. DATs CONREOUL Coordinate pairs (x,y). 


2.4 Program Structure 

Tiecelotd sr oOLtLanl program ConSisting of 11 subprograms 
interconnectei according to the abstract! aierarchical 
scheme in Figure A.2. The subprograns CONTOUR, SCAN, TRACE, 
CALC are modified versions of the nonIMSL routine called 
CONTUR, while the other ones were developed to the specific 
task describei in this text. 


2.5 Constraints and Assumptions 
The input file IMAG1.DAT is assumed to be an array of 
integer with dimension (256x256) with gray levels, i.e., 


pixel brightness value varying from zero to 255. 


1The reason to use the expression "abstract hierarchical 
scheme" is just to enforce hoteaethough Fortcan Language 
haS an intrinsic disjoint and nonnesting environmnent, 1.€., 
TettestiwiGnhilc, We Gan abstract an hierarchical structure 
for better understanding. 
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The faxinun number ™or Layers asec “= =2ze7 where tae 
Minus sign is used to indicate that the layers are to be 
chosen automatically by jividing the interval between the 
maximum and minimum pixel value into the absolute value of 
NC, rvather than with user's defined Layer values. 

The maximum number of points allowed in each countour is 
4800. This figure has been demonstrated more than sufficient 
during our experiment with aerialphotograph picture; never- 
theless, if any contour exceeds this limit the program will 
stop the contour in that point and an error message will be 
issued. The program will jump to the next contour, ‘and 
continue the 2xecution as if nothing have happened. 

This problem can be solved by just increasisg the dimen- 


Sion .Of X,Y NeGrons. 


2.6 Common Areas 
This program haS many common areaS Shared ody different 
subprograms. The table IX illustrates the relation between 


the subprogramns through these common areas. 


2.7 Error Messages 
Embedded inthe body of the respective coie listed in 


Appendix C. 


2.8 Usage 
The 11 subprograms that compose this program are grouped 
in a single module called THE8M.FOR. Thus, tiese are the 
following procsedupeS tO WOrKkK Wath theseplece cCfmeodc- 
Compilation: $FOR THES8M 
Linking 7 SLENKO? Bea 
Running : $RUN IHES8s 
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Pee GOgeanm Namwe; THESKNT.FOR 


EewObsgective 

Given a set of coordinates (x,/7) of the contour and the 
character matrix map yenerated by the program THESM.FOR, 
this program minimizes the number of knots in 2ach contour 
by selecting automatically the break points according to the 
Tagrange's minimization approach, i1.ée.: 

- Minimize the second derivative (maximize the smooth- 
ness) under the constraint of the maximum closeness error 
measurement less than a factor "S" which is a function of 
the number of breakpoints. For more details see thapter III. 

- The three-dimensional eifect is carried out by setting 
the weighting function W(I) proportional to the density of 
eS Mmrougacetn the neighborhood of each contour foint (x,y). 
These measurements are done by overlaping the contours 
record over the character matrix map which are shown in 


Pig ime A <3 . 


Seaerenput Data 

bate FORO11. DAT; record containing the reference and 
incremental w2ights, WO and AW respectively. 

File FORO18. DAT: Shabhactecumiatnix. Mapmeontaining the 
(Connectivity four) contours, represented by the letters: 
A,B,-«e-. Each letter defines a specific layer and the back- 
ground is represented by the blank character. 

File FORO20.DAT: record containing the (x,y) coordinates 


of the contours. 


Seo Output Data 
Pilewn@rnogewiAl: Crecocd iontaining the (x,y) coordinates 
of the selected knots of each contour. 


3.4 Pro@gtzam st ruceume 

This Proj mam consits @et Subprograns organized 
according to three abstract levels of hierarchy shown in 
Figure A.3. At the first level there is a main program which 
handles the input/output functions. At the s2cond level 
there are four subprograms which perform tae following 
functions: 

- Subprogram DPARAM: selection of the knots for the 
opened contours. 

~ Subprogram DPERCLO: selection of the knots for the 
closed contours. 

- Subprogram TRIDI: introduction of the weighting factor 
to the contours points according to the density of contours 
in the neighborhood of that area being analysed. 

- Subprogram DERIVATIVE: evaluation of th2 parametric 
functions "x" and "y" in the selected knots accomplished by 
the former subprograms. 

The third level programs perform fundamentil functions 
which descriptions are imbeded in the body of the source 


listing in Appendix C. 


3.5 Constraints and Assumptions 

The maximum number of coordinate pairs allowed per each 
contour is 4300. The character aap Must be a two- 
dimensional array with 256 rows and colunas. 

The following parameters were set for this experiment 
and they are assumed to be the default unless an explicit 
change 1S made. 

- WINDOW = 5: It gives the span of the adlgorithn “thas 
measures the density in the neighborhood of each point 
belonging to the contour. 

- K = 3: It gives the degree of the B-spliie; in this 


case they are cubic splines. 


—emeeAR =sJs It gives the option for automatis parametri- 
zation. This means that the parameter to represent each 
curve (contour) is to be constructed according to equation 
3.32. The reason to choose this proceiure is to enforce the 
robustness of the overall process, i.e., keep it as general 
aS possible. 
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pass mode. 


QO: This switch sets the program toa single 


=SNEST = 1000: It gives the maximum numb2r of knots 
allowed to be selected per each contour, what has denmon- 
strated enough for our experiments. 

—aret = O.VMem It is the requested pelative accuracy for 
Ehemwsooe Of F(p) =S. 

—Sepeic = 205 It 1s the maximum number 93£ iterations 
allowed in the estimation of the root of F(p) = S. 


3.6 Common Areas 

There are two common areas. The first one is named CIDA 
and it is used to share data between the main program and 
the subroutine TRIDI. The second common area 15 named OPT1 
and it is mentioned in the subroutines DPARAM and DPERCLO, 


but it is not used in our application here. 


3.7 Error Messages 
Embedded in the body of the respective cole listed in 


Ropendix C. 


3.8 Usage 
The 13 subprograms of this set are gGEOuUpea™ inmeowroun 


modules in the following fashion: 


MODULE SUBPROGKAMS 

THE8 KNI THESHANT (main} 

eae TRILL 

PPE RCD DPERC LO, ~ PENG 

DPARAM DPARAM, BSPLN, COSSIN, ROIS Gee 


BACK, KNOT, DISCO, RADPION, Deeper 
Each module iS compiled individually. The linking and 
running are performed in the following way: 
SLINK TIHESKNT, TRIDI, DPERCLO, DPR AM 
5 RUN een ee nein 


4. Program Nane: FILTREX.FOR 


4.1 Objective 

This program reads the file which contains the selected 
knots of the contours ani eliminates the repeated points 
(knots) which appear due to the round off imposed in the 
program THE8KNT.FOR during the storage of the selected pairs 
(X,Y) - 

Note: This routine eliminates the repetition pairs only 
inside each contour set not among different contours. This 


task is accomplished in the next Step. 


4.2 Input Data 

File FOROIR. DAT: (usually ER = 22) with {x,¥) coordi. 
nates of the selected knots (which may be repeated due to 
round off). The reason for truncation is that the resolution 


of the grid 15 17256"am a zero to 255 anteger range. 


5 Olepwt Data 
Ba ieee FORO DH. DAT: (asjaliy IW = 23) Melon) CoOoOrdi-— 
Mates of the selected knots without repetition in the same 


GQn0 ©COUr. 


424 Program Structure 
Only the Main program with no Subroutines or function 
Gael. 


4.5 Constraints and Assumptions 
The input records of selected knots are allowed to have 
at maximum 1000 elements. In order t> change this constraint 


it is necessary to change the dimension statement. 


4.6 Common Ar2as 


Not applicable in this case. 


4.7 Error Messages 
Embedded in the body of the respective cole listed in 


Appendix C. 


4.8 Usage 
Comeulation: SFOR FILT REX 
Lining >: $LINK FILTREX 
RUnMAng SoU) ie i Tags X 


5. Program Nate? RAStewee 


5.1 Objective 

This program reccnstructs the image as an i1niform grid 
of pixels with the same size as the original and presumably 
with as lower distortion as possible. The set 29f the scat- 
tered points (compressed data) is fed to an algorithm based 
on the triangulation that was discussed in Chapter IV. 


5.2 Lp ue vag 

File FORO19.DAT: number of contours in eacno level and 
the gray level value of that respective level. 

File FORO2Zs 20412 (77) coordinates of tae selected 


Knots. 


Ses OCULDUG Data 
File FORO31.DAT: unifocm (256x256) grid representing the 


reconstructed image ficture. 


IJ) ePrEOgram St pueeure 

This program consists of seven subprograms distributed 
in four layers. The Figura A.& Shows the abstract greranea. 
ical diagram. The description of each subprogram is provided 


ona 


in the source listing in the Appendix ¢c. 


5.5 Constraints and Assumptions 
The maxinum allowed number of input coordinate pairs 
(x,y) is 10000. 


5.6 Common Areas 
There is only one common area namely IDL?C, Which @ite 
Shared by the subprograms PWLOT and [DLCTN. 
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5.7 Error Messages 
Embedded in the body of the respective coie listed in 


Appendix C. 


pao Usage 
The seven subprograms are grouped into two modules which 


are organized in the following manner: 


MODULE SUBPROGRAM 
PRASM FRA5M 
PWLM1 EMGOL, SLD TANG, LDLCIN, 


PMLIN, EXIR, LoOXCHG 
Each module is compiled independently and th2 next steps 
are: 
Ponigias pot NK FRASM, PLM 
Runoing: SRUN FRASM 


6. Program Name; FRA4M.FOR 


6.1 Objective 

Analogous to the former program FRASM. FOR, this progran 
also reconstructs the uniform grid from a set of scattered 
points (selected knots). The difference here is regarding to 
the algorithm used to accomplish this task. Now a thin plate 
spline method is used to generate a C! surface, opposed to 
the former one that generates a C° surface. The descrirftion 


of this algorithm is in Chapter IV. 


Sezeeinput Data 

PasteerORQMZ. DAT: SparaneteretNPPR" that is the short for 
humber of points per region. This parameter sontrols the 
mesh of the auxiliary Jjcid over wnoich the thin plates 


splines are applied in order to reconstruct the image. 


a 


File FORO19. DAT: Nnumoer of COntOUrS Paeeds te ey eae 
the gray level value of the respective level. 
File FORO23.DAT: coorijiinate pairs (x,y) of the selected 


knots or pixels. 


6.3 Output Data 
File FORO31. DAT: Whe orm “(Zecwecy grid cepresenting 


reconstructed image picture. 


6.4 Program Structure 

This program consists of nine subprograns odrganized in 
five levels. The Figure A.5 shows the abstract hierarchical 
structure of the program tree. The description of each 


subprogram can be found in Appendix C. 


6.5 Constraints and Assumptions 
The maximum number of input coordinate pairs (xX, ¥). 2s 
15000. The moie switch was set to one, MODE = 1, that means: 
compute the coefficients for the local approximations by 
thin plate spline and return the grid of interpolated func- 
tion values nx0, x0, ny0O, yO?m@and eGmwhere;: 
nxO = 256 : number of rows 
x0 : row coordinate 
nyO = 256 : number of columns 
yO : column coordinate 


£0 >: gray level value of coordinat2 pair (x0,y9) 


6.6 Common Areas 


Not applicable in this caSe. 
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CopmeeEGOr Message 


bm@bedded in the body of the respective coi2 listed in 
Appendix C. 


6.8 Usage 
The nine subprograms are grouped into two addules which 


are organized in the following manner: 


MODULE SUBPROGRAM 
FRA4M FRA4M 
F24 LOPTS, EVLPTS, GRID 


POGCLer eC LOLe Ope V GRIT Ay 
DECCNC , woe LYE 
Each module is compiled independently and th® next steps 
ares 
Linking: SLINK FRAGH, F24 
Running; S$RUN FRA4M 


7. Program Name: DISP3M.FOR 


7.1 Objective 

This program reads tne (256x256) uniform grid (image) 
Progmunzt "LR™ (usually [TR = 31) and converts it intoa 
Gepeee S SOEmat picture, l.e., a 512x512 grid of pixels. In 


another words, an expanSion two by two is carriei out. 


7.2 Impwt Data 
Peer Onoin. DAL. = (USWai1,y —©e — 31) consisting of the 
uniform square grid (256x256) representing the image ina 


gray level range from zero to 256. 
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Je3 OEP see aed 
Pile RITA. DADsewecnststing of 312 reeords ara by te 


record suitable to be displayed in COMTAL. 


7.4 Porgtam SerwvGcedre 
Only the main program with no subroutines or functions 
Ca Jie 


7.5 Constraints and Assumptions 
AS we are dealing with pixels from the range zero to 255 
and due to the fact that Fortran uses a two's complement 
forn, each byte in Fortran Covers the gmterval ~ i\2eecomi- 
Therefore, in order to represent a range from zero to 255 
some transformation is accomplished in the following way: 
let "N" be the gray level for a particular pixel and "MB" be 
the actual element to be sent to the COMTAL device. The 
following algorithm will perform the necessary mapping: 
IF (N < 127) THEN 
MB = N 
ELSE 
MB 
END IF 


This procedure will guarantee that each pix21 value "N" 


-256 + N 


will be seen by the COMTAL aS an integer ranging from zero 
to 8252. 


7.6 Common Areas 
Not applicable in this case. 


7.7 Error Messages 
Embedded in the body of the respective code listed in 


Appendix C. 
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7.8 Usage 
Gemeitatzon ¢ SFOR DIS P3M 
Pan k 1 pg : SLINK DIS P3M 
Running : SRUN DIS P3M 


8. Program Name: RMSER 


8.1 Objective 
This projram gives a quantitative measuranent of the 
error between the original image (smooth version) and the 
reconstructed version from the compressed data. The metrics 
used are: 
- Root mean sguare error CRS) 
oer 1a 2 (FE, | - Fi )e/N2) 


el eo 


- Relative error ; co 
By ese oe (a Fi, )?/L 2 Fy, i 


es 


Geez Input Data 

Pate iIMAS1. DAT: consisting of the orljinal image 
(256x256) grid of gray level value in the range zero to 255 
to be used as reference. 

File FOROIR. DAT: (usually IR = 31) consisting of the 
uniform square grid (256x 256) of the reconstrusted image to 


be evaluated. 


Ses Outpet Data 
SePoote means Sduare CrrOrE  ; Eaysc 


- Relativ2 error : EReL 


8.4 Program Structure 
Only one block (main program). 


as 


8.5 Constraints and Assuaptions 

The input files IMAG1.DAT and FIRJ31.DAT ac2 assumed to 
be arrays of integer with dimensions (256x256) WLED  Gieay 
levels, i1.€., pixel brightness value varying €rom zero to 
Doo 


8.6 Common Areas 


Not applicable in this case. 


8.7 Error Messages 
Embedded in the body of the respective coie listed in 
Appendix C. 


8.8 Usage 
Compilatvonee. $F Chan eer 
Lamkang : SLU NKeeeeS Eek 
Running >: $RUN RMSER 
Nee” FLD eS 


As we can see from Figure A.1, the progranus and files 
are the kasic elements of our image procesing system. The 
files work iike a glue that links the programs among then- 
selves. In this section, the structure of th2 files and 
their contents are described. They 2re present2d in a chro- 
nological fashion as they are generated and used according 


to the diagram in Figure A.1. 


1s FLELE TADSGAM 
It is the input picture file Containing the  crrgriam 


image to be compressed. It is in the format suitable to be 


124 


displayed at COMTAL (image displaving device). The expan- 
Seon LactoOr LS fsur, Boot lhe 2 Gl td wwe eet > 12% oli2) 
Pueseteeoeeee this tile is the input of the program WOAT1.PAS and 
it is stored as a packed array. 


Peer tie IMNAG1.DAT 

This file is the output of the program WDAT1.PAS and the 
HIbSt Of them two inputs of pregram THESM.FOR. It contains 
image file in its reduced form, without the expansion factor 
four. Each pixel is represented by an "INTEGER*4+" type, in 
other words four bytes are reguired per each pixel. From now 
on, aS we will be dealing with the Fortran programs, the 
Fortran convension will be the default. For example, in the 
Petwmam Envibonment the generic file FORxyz.DAT is by 
default assigned to the logical unit xyz where x, y, andz 
ace digits from zero to nine. 

When we say that a file has the format 215, this means 
that each record of this file stores two integes> numbers in 
two consecutive fields of length five for each field, so 
that these numbers are right justified in theirc prespective 
fied dist 


Cea ee rOROT).DAT 

Memes tie SeCCOnd Input Of the program THESM-SFOR. This is 
i-camcloemput two parameters: "NC", number of coatour layers 
and "THRES", threshold value, so that the contours with less 
points than the threshold will be pruned. The format of the 
G@atarin thas file is 215. 


temo le FrOROQ?!].DAT 
Rietce al Cutout Of the program =~ THESKNT.F392 and it is 


used to input the two parameters WO, weight reference value 


5 


and AW, incremental weight. The format OF Vine 1a ea nee 
fil etwas Zone. 


Se, Ff ELE PORON De 

This file is an input of the program FRA44.FOR, thin 
plate algorithm implementation. It is used t) input the 
pacameter "NPPR" that is an estimate of the number of points 
pec region iato which the whole image will be divided and 
the thin plates evaluated. 

The fortatbmot this tikes ise. 


6. ELLE SOntibeeed 

This file 1s the character mapping accordinj to Figure 
3.4, and it contains the contour of the given image such 
that the layer is represented by a string of characters over 
a background nade up of blank characters. This file contains 
512 records, each one with the format 128A1. | 


Te Pfs FOR ONReD AT 

This file is the output of the program THE3IM. FORReeandg 
the input of the programs FOR4YM.FOIR and FORSM 1.FOR. Lt 
contains in e2ach record the information regarding to one 
Specific layer, namely the number of contours and the gray 
level value of such a layer. It will be used in the recon- 
stcuction prozsess. The foGmat of this £i2egs foal) ee 
sentinel record (-3, -3.) is used to mark the end of the 


data stream. 


8. FILE FOROZon Dar 
This file is an output of the program THESM.FOR and the 
input of the program THE8KNT. It contains the coordinate 


pairs (x,y) "Sf each “contour which are concatenated and 
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Separated by 2 tlag that gives the number of op2ints in the 
Zoerowrng GONEOUE and 1ts type (-?T 1f it is opened; -2 if it 
is closed.) 

Ghas Stats format 1S 275 and it is termimated by a 
Seatanel (=3,=3) . 


PoeeuLE POROZ2.DAT 

This file is the output of the program THESKNT.FOR and 
Pieme input of the program FILTREX.FOR. It sdntains the 
selected knots | (x,y) and its format is the same as 
Menoeue PA, 2.C., 215. 


ioe titer OROZS. OAT 

iimvoeciver!’S the output of the program FILITREX.DAT and 
the input of the prograns FOR4SM.FOR and FORTM.FOR. ie 
contains basically the same data as the file FOROQ22.DAT, 
except that the repeated pairs (x,y) in the same contour due 
tO the tEuncation from floating point to integec type repre- 
sentation are filtered out. 


iieeereeLe FOROS 1. DAT 

This file is the output of the program FRA4M.FOR and 
TeonerOoR, eadhad the input of the programs DISP3M.FOR and 
[open enw t CONtCalns the reconstpucted image, i.e., a grid 
O£ (256x256) integer values ranging from zero to 255, and 


ies eLOrMato1s 32T4. 


Ver eR TA. DA T 

Piersetilemisethe OouWtput Of the program DISP3M.FOR, and 
it 1s the reconstructed image suitable to be ijiisplayed at 
meme OMLAL , where some gualitative measurements are done. 
This file is already in the expanded form, i.e., (5.12512) 
Gieard % 


eZ 


In this Appendix A, we have seen the program environment 
and its two basic components: the programs ani the files. 
More information about the discussed programs are provided 
in the Appendix C where a source listing of each program can 


be found. 
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PROGRAM 
START 
INPUT FOROL0 
PICTURE OAT 
EIEE 


WOAT1 PAS FORO18 
DAT 


- THE 8M. FOR 








F ILTREX 
FOR 
{tH ) FORO12 FORO23 

DAT 
we ~—_. _. . ~~ FOR4EM.FOR 


pee eras ane) FORSM.F OR 





RMSER.FOR FORO3L 

(SO men an <8 ap os OO soe DAT 
OUNTITAT QUALLTAT DITA, DAT DISP3M.FOR 
MEASURE MEASURE 


( ENO 


Figure A.1 Program Environment. 
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Figure A.2 Contour Generating - Program Structure. 
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Figure A.3 
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Knot Selection - Program Structure. 
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IDXCHG 


Figure A.4 Triangulation Method - Program StEueture. 
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Figure A.5 Thin Plate Spline Method - Program Structure. 


Vers 


TASLE If 


Common Areas Cross References 







SUBPRO- 
GRAMS 
















CONTOUR 
COMP se 





TINA 


Bical 


ISABEL 


DAYHOF 


INTFAC 


TABS e e 
FEFE e ee «© e ee @ 
RITA e e 
Disks e @ e 


RIB 


a GEC ee ee ee eee ee 


134 


ee ee ee ee RN ee cee eens eel 


APPEND 
M 


MORE ON CO 


Ix 8B 
PRESSION 


The data compression scheme, as it was pres2nted in the 
former chapters, “is a two Step. procelure so thatatene first 
one is accomplished through the contour process and the 
second one oy the application of the knot selection algo- 
rithms over the contours generated in the first step. 
Besides, it is clear to see that this codinj process is 
irreversible. By this we mean that the original lata can not 
be exactly reconstructed from the compressed s2t of data, 
rather it can be only approximated, and as it was enforced 
in previous chapters, one of the goals of the whole image 
processing system was to reduce the intrinsic efcror associ- 
ated with the irreversible coding process. 

It was also seen that there is a trade-off between the 
achieved compression ratio and the error betwee1 the recon- 
structed imag? and the original one, 1i1.e., fdr a higher 
compression ratio we would expect higher error neasurements 
and valiee—wWemsi . 

There is another class of coding called reversible 
coding in which a set of given aumbers is napped into 
another set of numbers; the operation is error free and it 
fas Some important applications. A well known case of 
reversible code is the run length coding ijiscussed by 
Gonzales [Ref. 1]. 

A close examination of the compressed data jyenerated in 
the first and second phases shows that the stream of points 
representing the contour @2xhibits an interestiag property, 
the coordinate pair (x,y) and its successor or predecessor 


differs by a small number. 


S'S 


The Figure B.?t shows a histogram of the maximum incre- 
ments (x,y) for a typical contour diagram. AS we can See, 


the concentration of Stalyeraerementsets evicen 
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Figure B.1 Histogram for the Contours Increment. 


We can explore this characteristic by incceasing the 
compression ratio by a factor of about two in i reversible 
fashion with minimal effort. The process 1S somewhat similar 
to the run-length encoding and it works in the following 
way: insteai of using two bytes to store each selected 
knot, a (x,Y¥) pair of integer ranging from Zero to 255,55 
store only the first knot of each comtour ({origioj)  an@ fron 


there we store the increments between each point and its 


13.6 


predecessor. From the histogram in Figure 8.1, we can see 
that the majority of increments between consecitive points 
lies in the range from minus seven to plus seven, excluding 
Zero, because repeated points are not allowed in the 
contours (see function FILTREX in Appendix A). Taerefore, we 
can use the value correspondent to the increment zero asa 
flag to indicate that one of the two coordinates or eventu- 
ally both, have generated increment that has absolute value 
greater than seven. 

For each contour the selected knots are pasked in the 
EO leowrng way: 

~ First point (origin) uses two bytes: the First one to 
store x coordinate and the second one to store y coordinate. 

= From the second point and further on tae following 
rule is applied: if the absolute value of the increment of 
x and y are both less than seven, use one nible to store 
each ircrement; if one sf the coordinates or both have an 
absolute value greater than seven, than use the flag "zero" 
to indicate out of range followed by the value of increment 
of x stored in one byte and the increment of y stored in the 
next byte. 

The Figur2 B.2 shows the original stream of jata and the 
packed one. 

As we have seen, the algorithm is simple, efficient, 
easy to implenent, and it increases the data conpression by 
Geedetor Of about two without aincreaising the entropy of the 


System. 
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POR LFTEE*) ; 
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